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(4)  INTRODUCTION 

The  goal  of  the  proposed,  “Breast  Cancer  Screening  Using  Photonic  Technology”  research 
project  is  to  develop  optical  imaging  techniques  that  make  use  of  noninvasive  near-infrared 
(NIR)  light  for  obtaining  two-dimensional  (2-D)  transillumination,  and  three-dimensional  (3-D) 
tomographic  images  of  cancerous  lesions  of  human  breast.  The  imaging  method  involves 
illuminating  the  specimen  with  ultrashort  NIR  pulses  of  laser  light  and  construction  of  images 
using  two  approaches.  The  first,  known  as  the  shadowgram  method,  utilizes  the  image  bearing 
component  of  the  forward-transmitted  light  to  form  direct  shadow  images.  The  second,  referred 
to  as  the  three  dimensional  (3-D)  inverse  reconstruction  method,  makes  use  of  the  measured 
transmitted,  forward-scattered  or  backscattered  light  intensity  profiles,  known  or  estimated 
optical  properties  of  the  sample,  a  model  for  light  propagation  through  turbid  media  and  a 
sophisticated  computer  algorithm  to  construct  images  of  the  interior  structure  of  the  specimen. 

We  made  significant  advances  in  developing  both  of  these  approaches.  What  is  even  more 
important,  our  spectroscopic  imaging  experiments  with  normal  and  cancerous  human  breast 
tissues  reveal  useful  differences  in  optical  images  and  is  indicative  of  the  diagnostic  ability  of  the 
spectroscopic  imaging  method. 

(5)  BODY 

The  tasks  performed  and  the  progress  made  during  the  current  reporting  period  may  be 
broadly  grouped  as  follows: 

5.1  Enhancement  of  the  spectroscopic  imaging  arrangement, 

5.2  Time-sliced  and  spectroscopic  shadowgram  imaging  with  excised  human  breast  tissues 

5.3  Analytic  solution  of  Boltzmann  Transport  Equation,  and 

5.4  Development  of  the  inverse  image  reconstruction  method  using  backscattered  light. 

We  will  briefly  outline  our  accomplishments  in  each  of  these  areas,  and  refer  to  appended 

publications  (Appendices  1-4)  for  detailed  description  where  applicable. 

5.1  Enhancement  of  the  Spectroscopic  Imaging  Arrangement 

We  have  improved  on  the  spectroscopic  imaging  arrangement  that  we  assembled  and  used 
during  the  first  reporting  period  by  replacing  the  old  128x128  pixels  NIR  area  camera  with  a  new 
area  camera  equipped  with  a  320x240  pixels  sensing  element  (Sensors  Unlimited  SU-320).  The 
new  camera  has  improved  the  resolution  of  the  images  that  we  obtain  significantly.  This  is  an 
extension  of  the  Technical  Objectives  2  and  3  (TO  2, 3). 

5.2  Time-sliced  and  Spectroscopic  Shadowgram  Imaging  of  Excised  Human  Breast 
Tissues 

We  have  continued  to  pursue  the  time-sliced  and  spectroscopic  shadowgram  imaging  of 
excised  normal  and  cancerous  human  breast  tissues  using  the  experimental  arrangements 
developed  during  the  first  reporting  period  (TO  5-7,  Tasks  12, 14).  The  breast  tissue  specimens 
with  infiltrating  ductal  carcinoma,  and  infiltrating  lobular  carcinoma  from  patients  of  different 
ages  were  obtained  from  our  collaborators  at  the  Memorial  Sloan  Kettering  Cancer  Center  and 
National  Disease  Research  Interchange  under  an  IRB  approval  at  the  City  College  of  New  York. 

A  very  promising  and  interesting  result  of  2-D  spectroscopic  imaging  experiment  is  the 
wavelength-dependent  difference  in  light  transmission  through  the  cancerous  and  normal  tissues 
(TO  5-7,  Tasks  15).  As  a  measure  of  this  difference  we  monitored  the  ratio,  R  of  light  intensity 
transmitted  through  the  cancerous  tissue  to  that  through  the  corresponding  normal  tissue.  We 
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found  the  value  of  R  to  be  1.5  for  1225  nm  and  1.2  for  1300  nm,  a  significant  difference,  for  a 
breast  tissue  sample  with  poorly  differentiated  carcinoma,  grade  III,  with  sarcomatoid  features.1 
The  results  are  detailed  in  Appendix  1.  We  observed  similar  wavelength-dependent  variation  in  R 
for  breast  tissue  samples  with  ductal  carcinoma,  as  well.  We  will  pursue  measurements  involving 
normal  tissues  and  tissues  with  different  types  and  stages  of  cancer  to  examine  if  R  can  be  a 
parameter  whose  values  would  be  indicative  of  cancer. 

The  results  of  time-sliced  2-D  imaging  experiments  are  consistent  with  our  earlier  results2 
that  light  transited  through  the  cancerous  tissues  faster  than  through  the  normal  tissue. 
Consequently,  images  obtained  with  light  in  the  earlier  time  slices  highlighted  the  cancerous 
tissues,  while  those  obtained  with  light  in  the  later  time  slices  highlighted  the  normal  tissues. 
Time-sliced  imaging  can  thus  separate  out  normal  and  cancerous  tissues  in  excised  specimens. 
Some  of  these  results  are  also  presented  in  Appendix  1. 

5.3  Analytical  Solution  of  Boltmann  Transport  Equation 

Our  theoretical  endeavor  (TO  4,  Task  10;  TO  8,  task  18).  has  resulted  in  the  derivation  of 
analytical  solution  of  the  elastic  Boltzmann  Transport  Equation  in  an  infinite  uniform  isotropic 
medium  with  an  arbitrary  phase  function.3,4  This  new  approach  provides  a  more  accurate  and 
exact  description  of  photon  transport  through  highly  scattering  media  than  the  commonly  used 
diffusion  approximation  that  fails  to  adequately  account  for  ballistic  and  snake  photons.5  The 
approach  enables  calculation  of  (a)  the  exact  distribution  in  angle,  and  (b)  the  spatial  cumulants 
at  any  angle,  exact  up  to  an  arbitrary  high-order  ( Appendix  2)?  Terminating  the  cumulant 
expansion  at  the  second  order,  we  have  derived  an  analytical  solution  of  the  distribution  function, 
and  density  distribution  {Appendix  3).4  These  expressions  show  a  clear  picture  of  time  evolution 
of  particle  migration  from  ballistic  to  snake-like,  then  to  diffusive  regime.  Use  of  these  analytic 
solutions  will  provide  more  sophisticated  3-D  inverse  image  reconstruction  schemes. 

5.4  Development  Inverse  Image  Reconstruction  Method  using  Backscattered  Light 

We  have  made  another  major  advance  in  developing  a  novel  inverse  image  reconstruction 
(HR)  method  during  this  reporting  period  (TO  4,  Tasks  8,  9;  TO  8,  Task  18).  This  new 
approach6  uses  backscattered  photons  from  the  scattering  medium  containing  absorbing 
inhomogeneities,  and  differes  from  our  earlier  method7  that  used  transmitted  and  forward 
scattered  photons.  It  introduces  the  concept  of  propagation  of  spatial  Fourier  component  of  the 
scattered  wave  field  inside  the  scattering  medium,  and  then  develops  a  new  optical  diffuse 
imaging  methodology  based  on  that  theory.  The  method  is  fast  and  capable  of  providing  3-D 
image  information.  A  test  run  using  simulated  data  was  able  to  reconstruct  images  of  four 
inhomogeneities  located  up  to  2  cm  below  the  surface  of  a  human  tissue-like  semi-infinite 
scattering  medium  using  backscattered  photons.  The  inverse  reconstruction  method  using 
backscattered  light  is  detailed  in  the  preprint  of  a  paper  presented  here  as  Appendix  4.  We  will 
pursue  testing  of  this  method6  for  inverse  reconstruction  of  objects  using  experimental  data,  and 
compare  the  results  with  that  obtained  using  the  reconstruction  method  that  used  transmitted 
light. 7 

(6)  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Obtained  time-sliced  2-D  transillumination  images  normal  and  cancerous  tissues 
wherein  images  recorded  with  earlier  temporal  slices  of  transmitted  light  highlighted 
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cancerous  tissues  while  those  recorded  with  later  slices  accentuated  normal  fibrous 
tissues. 

•  Carried  out  spectroscopic  imaging  experiments  and  identified  a  ratio,  R  of  light 
intensity  transmitted  through  the  cancerous  tissue  to  that  through  the  corresponding 
normal  tissue  show  a  wavelength  dependent  variation  that  has  the  potential  to  be  used 
as  a  useful  parameter  for  cancer  identification. 

•  Developed  analytical  solutions  of  the  Boltzmann  transport  equation  that  enable  a 
more  accurate  description  of  the  ballistic  and  snake  components  of  light  emerging 
from  a  highly-scattering  medium  than  that  afforded  by  the  diffusion  approximation. 

•  Developed  the  theoretical  framework  and  computer  algorithm  for  inverse  image 
reconstruction  scheme  that  would  use  backscattered  light  to  provides  fast,  noise- 
resistant  3-D  images  of  objects  at  various  depths  inside  a  scattering  medium. 
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2.  W.  Cai,  M.  Lax,  and  R.  R.  Alfano,  “Analytical  solution  of  the  elastic  Boltzmann  transport 
equation  in  an  infinite  uniform  medium  using  cumulant  expansion,”  J.  Phys.  Chem.  B  104, 
3996  (2000). 
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(8)  CONCLUSIONS 

The  work  carried  out  during  this  reporting  period  affirms  some  of  our  earlier  inferences  and 
leads  to  some  new  conclusions.  First,  time-sliced  2-D  transillumination  images  recorded  with 
earlier  temporal  slices  of  transmitted  light  highlighted  cancerous  tissues  while  those  recorded 
with  later  slices  accentuated  normal  fibrous  tissues.  Second,  results  of  spectroscopic  imaging 
experiments  lead  to  a  ratio,  R  of  light  intensity  transmitted  through  the  cancerous  tissue  to  that 
has  the  potential  to  be  used  as  a  useful  parameter  for  cancer  identification.  Third,  analytical 
solutions  of  the  Boltzmann  transport  equation  that  we  obtained  enable  a  more  accurate 
description  of  the  ballistic  and  snake  components  of  light  emerging  from  a  highly-scattering 
medium  than  that  afforded  by  the  diffusion  approximation.  Fourth,  the  theoretical  formalism 
and  computer  algorithm  for  inverse  image  reconstruction  scheme  using  backscattered  light 
shows  (with  simulated  data)  the  potential  to  provide  fast  3-D  images  of  objects  at  various  depths 
inside  a  scattering  medium. 

“So  What  Section” 

The  implication  of  the  first  conclusion  is  that  time-sliced  imaging  offers  the  possibility  of 
highlighting  cancerous  lesions  in  human  breast.  The  second  conclusion  points  to  the  diagnostic 
potential  of  optical  imaging,  that  of  being  able  to  diagnose  a  tumor  as  it  is  being  imaged.  X-ray 
mammography,  most  often  used  method,  cannot  diagnose  cancer.  The  third  and  fourth 
conclusions  together  present  the  possibility  of  developing  robust  3-D  inverse  image 
reconstruction  formalisms,  that  in  addition  to  being  applicable  for  optical  mammography,  will  be 
useful  for  imaging  objects  inside  scattering  media,  such  as,  cloud,  fog,  smoke,  and  murky  water. 
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Proceedings  of  the  Inter-Institute  Workshop  on  In  Vivo  Optical  Imaging  at  the  NIH, 

Edited  by  A.  Gandjebache,  Optical  Society  of  America  (to  be  published),  p.142 

Temporally  and  spectrally  resolved  optical  imaging  of  normal  and  cancerous 

human  breast  tissues 

S.  K.  Gayen,  M.  Alrubaiee,  M.  E.  Zevallos  and  R.  R.  Alfano 

Institute  for  Ultrafast  Spectroscopy  and  Lasers,  New  York  State  Center  for  Photonic  Materials  and 
Applications,  Departments  of  Physics  and  Electrical  Engineering,  The  City  College  of  the  City  University 
of  New  York,  138th  Street  at  Convent  Avenue,  New  York,  NY  10031 
gaven@scisun.sci.ccny.cuny.edu,  alfano@scisun.sci.ccny.cuny.edu 

Abstract:  Time-sliced  and  spectroscopic  imaging  arrangements  were  used  to  obtain  two- 
dimensional  (2-D)  transillumination  images  of  a  composite  in  vitro  human  breast  tissue  sample 
comprising  cancerous  and  normal  fibrous  tissues,  adipose  tissues,  and  a  lymph  node.  Time-sliced 
imaging  approach  used  800-nm,  approximately  130-fs  duration,  1  kHz  repetition-rate  pulses  from 
a  Ti:sapphire  laser  system  to  illuminate  the  sample,  and  a  gated  imaging  system  that  provided  a 
variable-position,  ~80  ps-duration  electronic  gate  to  record  time-sliced  2-D  images.  Images 
recorded  with  earlier  temporal  slices  (approximately,  first  100  ps)  of  the  transmitted  light 
highlighted  the  lymph  node  and  cancerous  tissues,  while  the  later  slices  (later  than  300  ps) 
accentuated  the  adipose  and  normal  tissues.  Spectroscopic  imaging  arrangement  made  use  of 
1225  -  1300  nm  light  from  a  chromium-doped  forsterite  laser  for  sample  illumination,  a  Fourier 
space  gate  and  a  polarization  gate  to  sort  out  a  fraction  of  the  image-bearing  photons,  and  an 
InGaAs  area  camera  for  recording  2-D  images.  Marked  enhancement  of  image  contrast  between 
the  adipose  tissue  and  other  tissues  in  the  specimen  was  observed  when  the  wavelength  of 
imaging  light  was  near  resonant  with  the  1203-nm  optical  absorption  resonance  of  the  adipose 
tissue.  Wavelength-dependent  differences  in  relative  light  transmission  through  the  normal  and 
cancerous  tissues  were  observed. 

OCIS  codes:  (170.3880)  Medical  and  biological  imaging;  (170.6920)  time-resolved  imaging; 
(290.7050)  scattering,  turbid  media;  (170.6510)  spectroscopy,  tissue  diagnostics;  (170.3660)  light 
propagation  in  tissues;  (999.999)  Optical  mammography;  (999.999)  near-infrared  absorption 
spectroscopy  of  tissues;  (999.999)  spectroscopic  imaging;  (999.999)  time-sliced  imaging. 


Introduction 

Optical  mammography,  imaging  of  the  interior  structure  of  human  breast  with  light,  is  an  active  area  of 
optical  biomedical  imaging  research. [1-4]  Development  of  optical  breast  imaging  modalities  is  of  interest 
for  several  reasons.  Optical  imaging  methods  are  noninvasive  as  no  ionizing  radiation  is  involved.  Use  of 
different  wavelengths  of  light  has  the  potential  to  provide  diagnostic  information.  In  contrast  with  x-ray 
mammography,  light-based  methods  are  as  apt  to  image  dense  breast  of  a  younger  patient  as  that  of  an 
older  patient.  What  is  even  more  important,  inverse  image  reconstruction  methods  using  time-resolved  or 
frequency-domain  optical  measurements  may  provide  three-dimensional  (3-D)  tomographic  breast 
images.[5-7]  The  ability  to  generate  ultrashort  pulses  and  color  are  two  major  attributes  of  light  that  may 
be  exploited  to  develop  an  imaging  modality  with  diagnostic  ability. 

In  this  article,  we  present  the  results  of  time-sliced[7]  and  spectroscopic[8]  2-D  transillumination 
imaging  measurements  on  excised  human  female  breast  tissue  specimens  comprising  normal  and 
cancerous  tissues.  Time-sliced  imaging  makes  use  of  different  temporal  slices  of  the  transmitted  light 
to  form  2-D  images  following  the  illumination  of  the  sample  with  ultrashort  near-infrared  (NIR)  pulses 
of  light.  The  thrust  of  the  spectroscopic  imaging  experiment  is  to  examine  if  a  spectroscopic  difference 


may  be  used  to  enhance  image  contrast,  distinguish  between  different  types  of  tissues  in  a  specimen,  and 
obtain  diagnostic  information. 

Methods  and  Materials 

The  time-sliced  imaging  arrangement  used  800-nm,  approximately  130-fs  duration,  1  kHz  repetition-rate 
pulses  from  a  Ti:sapphire  laser  and  amplifier  system[9]  for  sample  illumination,  and  an  ultrafast  gated 
intensified  camera  system  (UGICS)  for  recording  two-dimensional  images  using  picosecond-duration 
slices  of  light  transmitted  through  the  sample.  The  UGICS  comprised  a  compact  time-gated  image 
intensifier  unit  fiber-optical ly  coupled  to  a  charge-coupled  device  (CCD)  camera.  It  provided  a 
minimum  gate  width  of  approximately  80  ps  whose  temporal  position  could  be  varied  in  25-ps  steps 
over  a  20-ns  range.  The  average  beam  power  used  in  the  experiment  was  approximately  200  mW.  The 
beam  was  expanded  by  a  beam  expander,  and  a  3-cm  diameter  central  part  of  it  was  selected  out  using 
an  aperture  to  illuminate  the  sample.  The  time-sliced  image  was  recorded  by  the  CCD  camera  and 
displayed  on  a  computer. 

The  experimental  arrangement  for  near-infrared  (NIR)  spectroscopic  imaging  made  use  of  1210- 
1300  nm  continuous-wave  mode-locked  output  of  a  Ct4+  :forsterite  laser  to  illuminate  the  sample.  A 
Fourier  space  gate[10]  in  conjunction  with  a  polarization  gate[ll]  selected  out  a  fraction  of  the  less- 
scattered  image-bearing  photons  from  the  strong  background  of  the  image-blurring  diffuse  photons.  A 
50  mm  focal-length  camera  lens  placed  on  the  optical  axis  at  a  distance  of  50  mm  from  the  aperture  in 
the  Fourier  gate  collected  and  collimated  the  low-spatial-frequency  light  filtered  by  the  aperture  and 
directed  it  to  the  128x128  pixels  sensing  element  of  an  InGaAs  NIR  area  camera.  The  average  optical 
power  of  the  incident  beam  was  maintained  at  approximately  35  mW  for  all  the  wavelengths  used  in  the 
imaging  experiment  using  appropriate  neutral  density  filters.  The  laser  beam  was  linearly  polarized 
along  the  horizontal  direction. 

The  composite  excised  breast  tissue  sample  used  in  the  experiments  reported  in  this  article  was 
assembled  from  tissues  obtained  following  the  modified  radical  mastectomy  of  a  30-year-old  patient.  It 
comprised  a  lymph  node  (LN)  with  surrounding  tissues,  a  piece  of  adipose  (A)  tissue,  and  a  piece  with 
normal  (N)  and  cancerous  (C)  fibrous  tissue.  Each  of  the  pieces  was  approximately  5  mm  thick,  and 
was  pressed  into  a  5-mm  thick  quartz  cell  to  ensure  uniform  sample  thickness  and  good  optical  contact 
between  the  adjacent  pieces.  According  to  an  accompanying  surgical  pathology  report,  the  cancer  was  a 
poorly  differentiated  carcinoma,  grade  III  with  sarcomatoid  features.  Figure  1(a)  shows  a  photograph  of 
the  exit  face  (the  side  that  faces  the  camera  in  the  experimental  arrangements  mentioned  above)  of  the 
sample  wherein  the  locations  of  different  types  of  tissues  in  the  composite  sample  are  tentatively 
labeled.  The  tissues  were  made  available  to  us  by  National  Disease  Research  Interchange  under  an  IRB 
approval  from  the  City  College  of  New  York. 

Results 

Time-sliced  Imaging 

Time-sliced  transillumination  images  of  the  sample  for  gate  positions  of  100  ps  and  350  ps  are  displayed 
in  Figs.  1(b)  and  1(c),  respectively.  The  zero  position  was  taken  to  be  the  time  of  arrival  of  the  light 
pulse  through  a  5-mm  thick  quartz  cell  filled  with  water.  The  spatial  intensity  profiles  of  the  images  in 
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Figure  l.(a)  A  photograph  of  the  exit  face  (the  side  that  faces  the  camera  in  the  experimental  arrangements)  of  the  composite 
breast  tissue  sample.  LN,  lymph  node;  A,  adipose  tissue;  N,  normal  fibrous  tissue;  C,  cancerous  tissue.  Time-sliced 
transillumination  images  of  the  sample  for  gate  delays  of  (b)  100  ps,  and  (c)  350  ps.  Spatial  profile  of  the  integrated 
intensity  distribution  along  a  horizontal  area  of  6  pixel  vertical  width  around  the  dashed  white  line  that  passes  through  the 
normal  fibrous  tissue  (dashed  line),  or  the  cancerous  tissue  (solid  line)  for  a  gate  delay  of  (d)  100  ps,  and  (e)  350  ps. 

Fig.  1(b)  and  Fig.  1(c)  integrated  over  two  6-pixel  wide  horizontal  areas  around  the  white  dashed  lines 
are  presented  in  Figs.  1(d)  and  1(e),  respectively.  The  two  areas  were  chosen  such  that  one  included  the 
normal  fibrous  tissue  in  the  upper  right  part  of  the  sample  while  the  other  included  the  cancerous  tissue 
in  the  lower  right  part  to  enable  close  comparison.  The  time-sliced  100-ps  image  clearly  highlights  the 
lymph  node,  adipose,  normal  fibrous,  and  cancerous  tissue  regions.  The  contrast  is  the  highest  between 
the  lymph  node  that  appears  the  brightest  and  the  adipose  tissue  that  appears  dark  in  the  image.  The 
spatial  intensity  distributions  of  the  100-ps  image,  displayed  in  Fig.  1(d),  show  the  highest  peak  in 
intensity  values  in  the  lymph-node  region  and  a  marked  trough  in  the  adipose  tissue  region  indicating 
much  higher  light  transmission  through  the  lymph  node  and  much  lower  transmission  through  the 
adipose  tissue  region  at  early  time.  More  interesting  is  the  contrast  between  the  cancerous  and  normal 
tissues  in  the  1 00-ps  image.  As  seen  in  the  right  side  of  the  image  and  the  spatial  intensity  profiles  of 
Fig.  1(d),  light  transmission  through  the  cancerous  tissue  is  significantly  higher  than  that  through  the 
normal  tissue. 

A  markedly  different  situation  is  observed  in  the  350-ps  image  of  Fig.  1(c)  and  the  corresponding 
spatial  intensity  profiles  of  Fig.  1(e).  The  adipose  tissue  region  appears  the  brightest,  and  the  spatial 
intensity  profile  peaks  in  the  adipose  tissue  region  indicating  much  higher  light  transmission  through  the 
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adipose  tissue  compared  to  transmission  through  other  tissues  in  the  sample  at  this  later  time.  What  is 
even  more  noteworthy,  the  difference  in  light  transmission  through  the  normal  and  cancerous  regions 
that  appeared  so  prominent  in  the  profiles  of  Fig.  1(d)  is  not  appreciable  at  this  late  time.  It  is  reflected 
by  the  close  overlapping  of  the  two  profiles  in  the  regions  of  the  normal  and  cancerous  tissues  in  the 
profiles  of  Fig.  1(e).  At  intermediate  times  (not  shown  in  figures)  relative  light  transmission  through  the 
lymph  node  and  cancerous  tissues  decreased  while  that  through  adipose  and  normal  fibrous  tissues 
increased  with  time.  Summarizing  the  time-dependent  transit  of  light,  we  find  that  the  light  transits 
fastest  through  the  lymph  node,  followed  by  that  through  cancerous  fibrous  tissue,  normal  fibrous  tissue, 
and  the  adipose  tissue.  Lower  scattering  or/and  higher  absorption  of  light  by  the  lymph  node  and 
cancerous  tissues  may  account  for  the  observed  temporal  behavior.  Since  there  is  no  known  absorption 
of  800-nm  light  by  breast  tissues,  we  attribute  these  time-dependent  differences  in  the  relative  light 
transmission  through  different  types  of  human  breast  tissues  to  the  differences  in  the  light  scattering 
characteristics  of  these  tissues. 

These  results  demonstrate  that  time-sliced  imaging  can  highlight  different  types  of  tissues  in  a 
sample.  What  is  even  more  important,  it  can  highlight  normal  fibrous  tissues  from  poorly  differentiated 
carcinoma  (grade  III)  with  sarcomatoid  features.  It  should  be  noted  that  more  pronounced  difference  in 
light  scattering  characteristics  between  normal  fibrous  tissues  and  infiltrating  ductal  carcinoma  was 
observed[12]  than  that  observed  in  this  study  between  normal  fibrous  tissues  and  poorly  differentiated 
carcinoma  (grade  III)  with  sarcomatoid  features.  This  in  turn  indicates  that  light  transport  characteristics 
will  vary  between  tissues  with  different  types  of  carcinoma,  as  well  as  between  normal  and  cancerous 
tissues. 

Spectroscopic  Imaging 

A  spectroscopic  difference  between  different  types  of  tissues  in  a  specimen  is  expected  to  provide  some 
distinguishable  signature  in  a  transillumination  image.  In  order  to  test  if  this  signature  may  be  realized  in 
practice,  we  obtained  images  of  the  sample  with  1225-nm  light  that  is  near-resonant  with  the  adipose 
tissue  optical  absorption  resonance  around  1203  nm,[13]  as  well  as,  with  light  of  wavelengths  away 
from  the  resonance.  Figures  2(a)  and  2(b)  show  a  'near-resonant  image'  recorded  with  1225-nm  light, 
and  a  typical  'nonresonant  image'  recorded  with  1300-nm  light,  respectively.  Figures  2(c)  and  2(d) 
display  the  corresponding  spatial  intensity  profiles.  The  solid  line  in  each  of  these  figures  shows  the 
profile  integrated  over  a  6-pixel  wide  area  around  the  long  dashed  line  that  runs  the  entire  length  of  the 
corresponding  image  and  includes  the  cancerous  tissue  region  in  the  lower  right  part  of  the  image.  The 
dashed  curve  superimposed  on  the  solid  curve  shows  the  profile  integrated  over  a  6-pixel  wide  t  normal 
fibrous  tissue  area  enclosed  by  the  smaller  box  in  the  corresponding  image.  The  solid  and  the  dashed 
curves  in  the  right  side  of  the  profiles  thus  enable  comparison  of  light  transport  characteristics  through 
normal  and  cancerous  tissues  in  the  specimen. 

The  salient  features  of  the  spectroscopic  images  and  corresponding  profiles  are:  (a)  the  adipose 
tissues  appear  much  darker  (less  light  transmission)  than  other  tissues  in  the  near-resonant  1225-nm 
image  as  compared  to  that  in  the  off-resonance  1300-nm  image;  (b)  cancerous  tissues  appear  brighter 
(higher  light  transmission)  than  the  normal  tissues  in  both  the  images;  (c)  while  the  overall  light 
transmission  through  the  normal  region  remains  approximately  at  the  same  level,  that  through  the 
cancerous  region  is  significantly  higher  at  1225  nm  than  at  1300  nm;  (d)  transmission  through  the  lymph 
node  exhibits  a  wavelength-dependent  variation  as  well,  being  higher  at  1225  nm  than  at  1300  nm.  The 
change  in  image  contrast  was  the  maximum  for  the  adipose  tissue  as  the  wavelength  of  imaging  light 
was  changed  from  1225  to  1300  nm.  Adipose  tissue  region  appeared  as  a  much  deeper  trough  in  the 
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Figure  2.  Spectroscopic  2-D  transillumination  image  of  the  breast  tissue  sample  described  in  the  text  for  light  of  wavelength  (a)  1225  nm, 
and  (b)  1300  nm.  Corresponding  spatial  profiles  are  shown  in  (c)  and  (d),  respectively.  Spatial  profiles  are  the  integrated  intensity 
distribution  along  a  6-pixel  wide  horizontal  area  around  the  white  dashed  line  covering  the  entire  length  of  the  sample  including  the 
cancerous  tissue  region  (solid  line  in  the  profiles)  and  that  around  the  small  dashed  line  denoting  the  normal  tissue  region  (broken  line  in 
the  profiles). 


spatial  intensity  profile  of  the  1225-nm  image  compared  to  that  for  the  1300-nm  image.  For  a  more 
quantitative  description  of  the  observed  behavior,  we  monitored  the  image  contrast,  C(X)  =  {IF  -IA)/(If 
+IA),  where  IA(X)  is  the  optimal  intensity  value  at  wavelength  X  on  the  spatial  profile  of  the  image  at  the 
adipose  tissue  location,  and  Ip{X)  is  the  corresponding  intensity  in  the  immediate  fibrous  tissue  region. 
Value  of  contrast  at  1225  nm  is  0.27  and  0.10  at  1300  nm.  As  the  laser  output  was  tuned  away  from 
1225  nm  to  off-resonance  wavelengths,  the  contrast  between  the  adipose  and  fibrous  regions  in  the 
images  decreased  from  that  maximum  value  of  0.27  towards  0.10.  These  results  clearly  demonstrate 
that  an  appreciable  spectroscopic  difference  may  significantly  enhance  the  contrast  between  different 
types  of  breast  tissues  in  a  transillumination  image,  and  is  consistent  with  our  previous  results  with 
adipose  and  fibrous  human  breast  tissues.  [8] 

Even  more  promising  and  interesting  is  the  wavelength-dependent  difference  in  light 
transmission  through  the  cancerous  and  normal  tissues.  As  a  measure  of  this  difference  we  may  monitor 
the  ratio,  R  of  light  intensity  transmitted  through  the  cancerous  tissue  to  that  through  the  corresponding 
normal  tissue.  Taking  the  averaged  intensity  values[14]  around  the  middle  of  the  normal  and  cancerous 
tissues  (Pixel  #  1 10  in  figures  2(c)  and  2(d)),  we  obtain  the  value  of  R  to  be  1.5  for  1225  nm  and  1.2  for 
1300  nm,  a  significant  difference.  We  observed  similar  wavelength-dependent  variation  in  R  for  ductal 
carcinoma  and  normal  breast  tissue  samples  as  well.  More  measurements  involving  normal  tissues  and 
tissues  with  different  types  and  stages  of  cancer  are  needed  to  examine  if  R  can  be  a  parameter  whose 
values  would  be  indicative  of  cancer. 

In  summary,  the  spectroscopic  and  time-sliced  imaging  methods  show  tissue  selectivity.  A 
combined  spectroscopic  and  time-sliced  imaging  approach  has  the  potential  to  provide  more  information 
even  than  the  x-ray  techniques. 
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We  study  the  analytical  solution  of  the  time-dependent  elastic  Boltzmann  transport  equation  in  an  infinite 
uniform  isotropic  medium  with  an  arbitrary  phase  function.  We  calculate  (1)  the  exact  distribution  in  angle, 
(2)  the  spatial  cumulants  at  any  angle,  exact  up  to  an  arbitrary  high  order  n.  At  the  second  order,  n~  2,  an 
analytical,  hence  extremely  useful  combined  distribution  in  position  and  angle,  is  obtained  as  a  function  of 
time.  This  distribution  is  Gaussian  in  position,  but  not  in  angle.  The  average  center  and  spread  of  the  half¬ 
width  are  exact.  By  the  central  limit  theorem  the  complete  distribution  approaches  this  Gaussian  distribution 
as  the  number  of  collisions  (or  time)  increases.  The  center  of  this  distribution  advances  in  time,  and  an  ellipsoidal 
contour  that  grows  and  changes  shape  provides  a  clear  picture  of  the  time  evolution  of  the  particle  migration 
from  near  ballistic,  through  snake-like,  and  into  the  final  diffusive  regime.  This  second-order  cumulant 
approximation  also  provides  the  correct  ballistic  limit.  Algebraic  expressions  for  the  nth  order  cumulants  are 
provided.  The  number  of  terms  grows  rapidly  with  n ,  but  our  expressives  are  recursive  and  easily  automated. 


I.  Introduction 

Search  for  an  analytical  solution  of  the  time -dependent  elastic 
Boltzmann  transport  equation  has  lasted  for  many  years.1-3 
Besides  being  considered  as  a  classical  problem  in  fundamental 
research  in  statistical  dynamics,  a  novel  approach  to  an  analytical 
solution  of  this  equation  may  have  applications  in  a  broad  variety 
of  fields.  To  our  knowledge,  an  exact  solution,  even  in  an  infinite 
uniform  medium,  is  available  only  for  isotropic  scattering  case, 
given  by  E.  H.  Hauge,4  in  the  form  of  a  Fourier  transform  in 
space  and  Laplace  transform  in  time.  Based  on  the  angular 
moment  expansion  with  cut-off  to  certain  order,  the  Boltzmann 
transport  equation  is  transferred  to  a  series  of  moment  equations. 
In  the  lowest  order,  a  diffusion  equation  is  derived  and  its 
analytical  solution  in  an  infinite  uniform  medium  is  obtained 
for  anisotropic  scattering  cases.  This  analytical  solution  has  been 
broadly  applied  in  many  applications.  For  example,  the  solution 
of  inverse  problems  in  optical  tomography,  such  as  the  location 
of  a  tumor  in  a  woman's  breast  from  the  scattering  of  light 
pulses,  requires  the  inversion  of  a  weight  matrix5  obtained  by 
convoluting  two  Green's  functions  of  the  forward  scattering 
problem.  The  analytical  solution  of  the  diffusion  equation  has 
provided  the  needed  Green's  function.  A  similar  procedure  can 
be  applied  to  other  problems,  such  as  using  a  laser  to  monitor 
cloud  distributions,  to  detect  objects  inside  a  cloud,  or  the  use 
of  low-frequency  sound  to  detect  oil-bearing  layers  deep  under 
water.  The  diffusion  approximation  fails  at  early  times  when 
the  particle  distribution  is  still  highly  anisotropic.  The  solutions 
of  the  diffusion  equation  or  the  telegrapher's  equation  do  not 
produce  the  correct  ballistic  limit  of  particle  propagation.6 
Numerical  approaches,  including  the  Monte  Carlo  method,  are 
the  main  tools  in  solving  the  elastic  Boltzmann  equation; 
however,  detailed  solution  of  a  five-dimensional  Boltzmann 
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transport  equation  using  a  predominately  numerical  approach 
leads  to  prohibitive  CPU  times. 

In  this  paper,  we  seek  an  analytical  solution  of  the  elastic 
Boltzmann  transport  equation  in  an  infinite  uniform  medium. 
We  assume  that  the  phase  function,  P( s,  so),  depends  only  on 
the  scattering  angle:  P( s,  so)  =  P(s-s0),  where  the  velocity  v  = 
us,  s  is  a  unit  vector  of  direction,  and  v  is  the  (constant)  speed 
in  the  medium.  Under  this  assumption,  we  can  handle  an 
arbitrary  phase  function.  We  obtain  the  exact  angular  distribution 
as  a  function  of  time.  Based  on  this  solution,  we  use  a  cumulant 
expansion  of  the  particle  distribution,  7(r,  s,  /),  and  derive  exact 
spatial  cumulants  up  to  an  arbitrary  high  order  at  any  angle 
and  time.  A  cut-off  at  second  order  yields  a  simple  analytical 
expressions  for  7(r,  s,  t),  as  a  function  of  position  r,  angle  s, 
and  time  /,  and  the  particle  density  distribution,  N( r,  t),  as  a 
function  of  position  r  and  time  t.  These  spatial  Gaussian 
distributions  have  the  exact  first  cumulant  (the  position  of  center 
of  the  distribution)  and  the  exact  second  cumulant  (the  half¬ 
width  of  spread  of  the  distribution).  After  many  scattering  events 
have  taken  place,  the  law  of  large  numbers  (the  central  limit 
theorem)  guarantees  that  the  spatial  Gaussian  distribution  that 
we  calculate  will  become  accurate  in  detail,  since  the  higher 
cumulants  become  relatively  small.  At  early  times,  the  spread 
of  the  distribution  is  narrow,  hence,  the  spatial  distribution 
function  can  be  claimed  quantitatively  accurate  for  many 
applications,  in  the  sense  that  it  has  the  correct  mean  position 
and  the  correct  half-width  of  spread  as  a  function  of  time. 
Measurement  of  the  higher  order  cumulants  could  require 
measuring  instruments  of  extreme  resolution. 

The  remainder  of  this  paper  is  organized  as  follows.  Section 
II  describes  the  derivation  of  formula:  (1)  obtaining  an  exact 
solution  of  the  distribution  in  angle,  (2)  obtaining  an  exact 
formal  solution  in  position  and  angle,  (3)  using  the  cumulant 
expansion  to  calculate  the  exact  analytical  expressions  of 
cumulants  up  to  an  arbitraiy  high  order,  (4)  describing  the 
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calculation  of  the  particle  distribution  function  using  a  spatial 
Fourier  transform.  Section  III  discusses  using  a  cut-off  at  second 
order  to  produce  explicit  expressions  of  the  distribution  function 
and  the  density  distribution.  A  brief  discussion  and  summary 
then  follows  in  Section  IV.  In  the  Appendix,  we  derive  analytical 
formulas  for  evaluating  integrals  in  eq  12. 

II.  Derivation  of  Cumulants  to  an  Arbitrary  High  Order 

The  clastic  Boltzmann  kinetic  equation  of  particles,  with 
magnitude  of  velocity  v ,  for  the  distribution  function  /( r,  s,  /) 
as  a  function  of  time  /,  position  r,  and  direction  s,  in  an  infinite 
uniform  medium,  from  a  point  pulse  light  source,  <$(r  —  ro)  <5(s 
—  so)  <5(f  —  0),  is  given  by3 

dl( r,  s,  t)/dt  +  vs*  V^r,  s,  t)  +  fij( r,  s,  t)  = 

fij  F(s,  s7)[/(r,  s',  t)  -  /( r,  s,  t)]ds'  + 

<5(r  —  r0)  <5(s  —  s0)  A(f  —  0)  (1) 

where  fts  is  the  scattering  rate,  jUa  is  the  absorption  rate,  and 
P(s ,  s)  is  the  phase  function,  normalized  to  /  ds7  P{ s',  s)  =  1. 
When  the  phase  function  depends  only  on  the  scattering  angle 
in  an  isotropic  medium,  we  can  expand  the  phase  function  in 
Legendre  polynomials  with  constant  coefficients, 


s,  /)  (the  source  is  located  at  ro  —  0)  is  given  by 

/( r,  s,  0  =  (<5(r  -  v  f‘0s(t'W)  <5(s(t)  -  s)>  (5) 

where  (...)  means  the  ensemble  average  in  the  velocity  space. 
The  first  6  function  imposes  that  the  displacement,  r  —  0,  is 
given  by  the  path  integral.  The  second  6  function  assures  the 
correct  final  value  of  direction.  Equation  5  is  an  exact  formal 
solution  of  eq  1 ,  but  can  not  be  evaluated  directly.  We  make  a 
Fourier  transform  for  the  first  ^-function  in  eq  5,  then  make  a 
cumulant  expansion,8  and  obtain 


/( r,  s,  t)  =  F( s,  s0,  t) - f  dk  exp]  /k*r  + 

(2ny  [ 


(2jz) 

n=l  nl  jn  )\ 


(6) 


where  T  denotes  time-ordered  multiplication.9  In  eq  6,  index  c 
denotes  cumulant,  which  is  defined  in  many  statistics  text¬ 
books,10  as  (A)c  =  (A),  (A2)c  —  (A2)  -  (A)(A),  and  a  general 
expression  relating  (Am)  and  (Am)c,  which  is  given  by: 


P(s,  s')  =  ~^LalPl( s’s')  (2) 

We  first  study  the  dynamics  of  the  distribution  in  direction 
space,  F( s,  so,  t),  on  a  spherical  surface  of  radius  1,  which  is 
equivalent  to  the  velocity  space  in  the  elastic  scattering  case. 
The  kinetic  equation  for  F( s,  So,  t)  can  be  obtained  by  integrating 
eq  1  over  the  whole  spatial  space,  r.  The  spatial  independence 
of  yUs,  and  P( s,  s')  retains  translation  invariance.  Thus  the 
integral  of  eq  1  obeys 


s±(«w^r... 

j!\  1!  I  m2!\  2!  /  mn\\  n\  I 
6(m  —  ml  ~  lm2  —  ...  nmn  —  ...)  (7) 

Hence,  if  (Am)  m  —  1,  2,  ...  n  have  been  calculated,  (Am)c  m  = 
1,  2,  ...  n  can  be  recursively  obtained  and  conversely.10  In  the 
following,  we  derive  the  analytical  expression  for  the  ensemble 
average  (f{Qdtn  ...  f‘0dti  T[sjn(tn)  ...  fy(*i)]>.  Using  a  standard 
time-dependent  Green’s  function  approach,  it  is  given  by 


9F(s,  s0,  t)/dt  +  jUaF(s,  s0,  t)  +  fis[F( s,  s0,  t)  - 

/ P(s,  s')  F( s',  s0,  t)ds']  =  <5(s  -  s0)  6(t  -  0)  (3) 

In  contrast  to  eq  1,  if  we  expand  F( s,  So,  t)  in  spherical 
harmonics,  its  components  do  not  couple  with  each  other. 
Therefore,  it  is  easy  to  obtain  the  exact  solution  of  eq  3:7 

^2 1+  1 

F( s,  s0,  t)  =  expC-^OX - cxpi—gfiP {(s*s0) 

l  4  71 

=  exp(-,Ma0X - exp(-ft/)^y,m(s)}^,(s0)  (4) 

/ 

where  gi  —  (is[  1  —  ail (21  +  1)].  Two  special  values  of  gi  are  go 
=  0,  which  follows  from  the  normalization  of  P( s,  s7)  and  g\  = 
v/lu  where  k  is  the  transport  mean  free  path,  defined  by  lt  = 
v![fis(  1  cos0)],  where  cos 6  is  the  average  of  s*s7  with  P( s,  s') 
as  weight.  In  eq  4,  F/m(s)  are  spherical  harmonics.  Equation  4 
serves  as  the  exact  Green’s  function  of  particle  propagation  in 
velocity  space.  Since  in  an  infinite  uniform  medium  this  function 
is  independent  of  the  source  position,  ro,  requirements  for  a 
Green’s  function  are  satisfied,  especially,  a  Chapman— Kol¬ 
mogorov  condition  is  obeyed:  /  ds'F(s",  s7,  t  —  t')  F( s7,  s,  /'  — 
f0)  =  F(s",  s,  t  -  t0).  In  fact,  in  an  infinite  uniform  medium, 
this  propagator  determines  all  particle  migration  behavior, 
including  its  spatial  distribution,  because  displacement  is  an 
integration  of  velocity  over  time.  The  distribution  function  7(r, 


J ds(1)F(s,  s(n),  /  -  tn)s<*Wn\  s(n-l\  tn  -  tn_1)4nr]1> ... 


F(sa\  s(1),  t2  -  ty;’F(s-\  So,  r,  -  0)  +  perm}  (8) 


.0)w.O) 


where  the  word  “perm”  means  all  n!  —  1  terms  obtained  by 
permutation  of  {/,•},  i  ==  1,  ...  n,  from  the  first  term.  In  eq  8, 
F( s(,),  s(,_1),  ti  —  t,-i)  is  given  by  eq  4.  Since  eq  4  is  exact,  eq 
8  provides  the  exact  nth  moments  of  the  distribution.  In 
Cartesian  coordinates  three  components  of  s  are  [s*,  sy ,  sz].  For 
convenience  in  calculation,  however,  we  will  use  the  compo¬ 
nents  of  s  on  the  basis  of  spherical  harmonics: 


s  —  [jj,  Sq ,  *s,_j]  =  [Tjj(s),  T10(s),  T|_j(s)]  — 

f— 2_1/2sin  6  e+i*,  cos  0,  +2~m  sin  6  e^] 

The  recurrence  relation  of  the  spherical  harmonics  is  given  by 

YJi* 1.  m,j\l  +  i,  m  +  j)  x 
/ 

</,  lf0,0|/  +  i,0>,i  =  ±l  (9) 

where  (fi,  h,  wi,  mi\l,  m)  is  the  Clebsch— Gordan  coefficients 
of  angular  momentum  theory,11  which  are 
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{I  —  i,\,mj\l,m  +  j)  = 


*(/  -  m)(l  -  m+  1)' 

1/2 

(/  +  m)(l  —  m  +  1)' 

1/2 

(l  +  m)(l+m  +  1)' 

(21-1)21 

21(1  -h  1 ) 

.  (21  +  2)(2l  +  3)  . 

'(/  -  m)(l  +  m) 

1/2 

m2  1 

1/2 

(l  +  m+  1)(l  —  m  +  1)' 

(21-1)1  \ 

W+  1)J 

(l  +  1  )(2/  +  3)  J 

'(/  +  m)(l  +  m  + 

1/2 

(/  —  m)(l  -F  m  4-  1) 

1/2 

’(/  —  m)(l  —  m  +  1) 

(2/- 1)2/ 

2/(/+l) 

.  (21  +  2)(2/  +  3)  . 

with  the  row  index  (from  above)  y  —  — 1,0,  1  and  the  column 
index  (from  left)  i  =  1,  0,  -1.  The  orthogonality  relation  of 
spherical  harmonics  is  given  by 

/  ds  >7v(s)  K/„/s)  =  27^-jA  ( 1  °) 

Using  cqs  9  and  10,  integrals  over  ds(n)  ...  ds(1)  in  eq  8  can  be 
analytically  performed.  We  obtain,  when  so  is  set  along  z,  that 


</o'd?n-/o'd',n^„).^|('1)]>  = 

2(/-Z'm)+1 

\  m—\  / 

-  X 

4;r 

n  n-k+ 1  *-l 

n<‘ -  x  >• 

n-A  k  n-k+ 1 

X  ^O.  01 

m=l  m— 1  1 


1 


F(s,  s0,  t) 


X%-  ,/jwX-X 


The  exact  second  cumulant  provides  the  correct  half-width  of 
spread  of  the  distribution.  Moreover,  the  central  limit  theorem 
claims  that  as  the  number  of  collision  events  become  large 
enough,  the  resulting  Gaussian  distribution  approaches  detailed 
accuracy  beyond  first  two  exact  cumulants.  At  early  time,  spread 
of  the  spatial  distribution  is  narrow,  possibly  narrower  than  the 
available  detection  instruments,  hence,  a  spatial  distribution  with 
exact  first  and  second  cumulants  may  provide  an  accurate 
enough  description  of  particle  distribution  for  many  applications. 

For  the  reader’s  convenience,  the  expressions  below  are  given 
in  Cartesian  coordinates  with  indices  a,  (3  =  [*,  y,  z]-  These 
expression  is  obtained  by  use  of  an  unitary  transform  -ya  —  Uaftj 
j  =  1,  0,  -1  from  eq  11,  (up  to  second  order)  which  is  based 
on  sj  ~  Ti/s),  with 


U  = 


—2~m  0 
2~mi  0 
0  1 


We  set  So  along  the  z  direction  and  denote  s  as  ( 0 ,  <p).  Our 
cumulant  approximation  to  the  particle  distribution  function  is 
given  by 


l-  X'»-m+l’°>|+Perm  (11) 

m— 1 


Hr,  s,  t)  = 

F( s,  s0,  0 


(4 7ifn  (del B)1/2 


exp|-~(fl  l)^(r  -rc)a(r  -  r^j  (13) 


where  i/=  ±1,/=  1,  2,  ...n,  and 

exp(-^fll){  fodtnfodtn-l  fodt'  CXP^S,(t  ~  01 

exP[  S/— i„(f«  “  _  °)]}  U2) 

Note  that  all  ensemble  averages  have  been  performed.  Equation 
12  involves  integrals  of  exponential  functions,  which  can  be 
analytically  performed.  An  explicit  expression  for  evaluating 
integrals  in  eq  12  is  presented  in  the  Appendix.  Equation  12 
includes  all  related  scattering  and  absorption  parameters,  gi ,  l 
=  0,  1,  ...  and  /xa,  and  determines  the  time  evolution  dynamics. 
The  final  particle  direction,  s,  appears  as  argument  of  the 
spherical  harmonics  y/m(s)  in  eq  11.  Substituting  eq  12  into  eq 
1 1,  and  using  a  standard  cumulant  procedure,  the  cumulants  as 
functions  of  angle  s  and  time  t  up  to  an  arbitrary  nth  order  can 
be  analytically  calculated.  The  final  position,  r,  appears  in  eq 
6,  and  its  component  can  be  expressed  as  |r|  Ki/f),  j  —  1,  0, 
—  l,  with  |r|  and  r  are,  separately,  magnitude  and  unit  direction 
vector  of  r.  Then,  performing  a  numerical  three-dimensional 
inverse  Fourier  transform  over  k,  an  approximate  distribution 
function,  7(r,  s,  f),  accurate  up  to  nth  cumulant,  is  obtained. 

IIL  Gaussian  Approximation  of  the  Distribution  Function 

By  a  cut-off  at  the  second  cumulant,  the  integral  over  k  in 
eq  6  can  be  analytically  performed,  which  directly  leads  to  a 
Gaussian  spatial  distribution  displayed  in  eq  13.  The  exact  first 
cumulant  provides  the  correct  center  position  of  the  distribution. 


with  the  center  of  the  packet  (the  first  cumulant),  denoted  by 
r^,  located  at 

r\  =  G^A,P,(cos  0)[(l  +  -  gl+])  +  l fig,  ~  ft.,)] 

(14.1) 

rcx  =  tiyjitfX cos  0)  cos  4>\figt  -  g,~  i)  -fig,-  ft+i)] 

(14.2) 

where  G  =  v  exp(-fiat)/F(s,  so,  t).  A,  =  (l/4;r)  exp(-£,r),  gi  is 
defined  after  eq  4,  and 

fig)  ~  [exp(gf)  —  l]/g  (15) 

ry  is  obtained  by  replacing  cos  <f>  in  eq  14.2  by  sin  (p .  In  eq  14, 
/>[m)( cos  6)  is  the  associated  Legendre  function. 

The  square  of  the  average  spread  width  (the  second  cumulant) 
is  given  by 


Bafi  —  vGA^  —  rcarCp/2  (16) 

where  all  the  coefficients  are  functions  of  angle  and  time: 


a*=Xa'p/(cos  6) 


"-y +(,+i)(,+y+ 


21-  1 


21-  1 


■£<3,+ 


21  +3 
(!+l): 


2Z  +  3 


•£<4> 


(17.1) 
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A«.yv  =  'L~AtPi(C0S  0) 

I  2 


/(/“!) 


2/- 


</+1*/+y + /(/- y + (/+1*/+y 

21  +  3  '  2/  -  1  '  2/  +  3  '  . 


y.-VfVeos  0)  cos(2  </>) 
/  2 


1  -£5"  +  — — fif’- 


2/  -  1  2/  +  3 


-£<3> - £?> 


(17.2) 


2/  -  1  2/  +  3 

where  (+)  corresponds  to  and  (— )  corresponds  lo  Avv 

1 


A„ =  Avi  -  Y-A,P{2,(cos  0)  sin(20) 
/  2 


2/-  1 


•£<1)  + 


1  -£® - — £<3> - — £}4> 


2/  +  3  2/  -  1 


,1 


2/ +  3 


(17.3) 


A,z  =  Az*  =  6)  cos  0; 

i  2 


2(/-  1) 


2/-  1 


-£<1)- 


—  +  2)E^2>  +  — — £f>  +  — — £j4) 

2/  +  3  2/  —  1  2/  +  3  . 


(17.4) 


A yz  is  obtained  by  replacing  cos  <p  in  eq  17.4  by  sin  (p .  In  eq 
17.1-17.4 


4' =m- ft- 2) -fa, - - 8,-2)  os.!) 

=\fig,~  8m)  -fig,  ~  8,+i)V(8,+\  ~  8m)  (18-2) 

^3)  =  [/(ft-ft-,)-*]/(ft-ft-,)  (18-3) 

^4)  =  Wft-ft+I)-f)/(S/-ft+ 1)  (18.4) 

A  cumulant  approximation  for  the  particle  density  distribution 
is  obtained  from  the  exact  expression  N( r,  t)  =  (d(r  —  vf!0s(t') 
d t')).  Using  /ds  F(s,  s',  t)  =  exp(—  fiat),  we  have  a  Gaussian 
shape 


N(r,  t) : 


(AjlD^Vtf2  A*DxxVt 


exp 


{z-RX 


4  Duvt 


exp 


(x  2  +  /) 


4Drrvt 


exP(— /V)  (19) 


with  a  moving  center  located  at 

Rz  =  v[l  “  exp(-g,/)]/g,  (20) 

and  the  corresponding  diffusion  coefficients  are  given  by 
^  v  \  t  3ft  _  82  rl  ,  , 

— — - :[1  -  exp(-g2/)]  -  2_[i  -  exp(-g,0)2[  (21.1) 

fttft  “  ft)  2g3 


— 11  -  exp(-g,/)l  - 

ft) 


[1 -exp(-g20)  (21.2) 

ft(S  i  ~  ft)  J 

In  contrast  to  eqs  14  and  17,  these  results  are  independent  of 
gi  for  /  >  2.  Figure  1  shows  the  moving  center  of  particles,  Rz 
(eq  20),  and  the  diffusion  coefficients,  D a  and  D ^  (eqs  21),  as 
a  function  of  time,  where  g/  is  calculated  by  Mic  theory  of  light 
scattering12  assuming  (for  this  figure)  that  the  “uniform” 
scattering  medium  consists  of  water  droplets  with  r/X  =  1  arc 
uniformly  distributed  in  air,  with  r  the  radius  of  the  droplet,  X 
the  wavelength  of  light,  and  the  index  of  refraction  m  =  1.33. 

Each  distribution  in  eqs  13  and  19  describes  a  particle  “cloud” 
anisotropically  spreading  from  a  moving  center,  with  time- 
dependent  diffusion  coefficients.  At  early  time  t  — *  0,  f(g)  %  / 
+  0(t2)  in  eq  15,  and  Efp  ~  t2l  2  +  0{t2)  for  j  —  1,  2,  3,  4  in  eqs 
18.  From  eqs  14,  17,  and  20—21,  we  see  that  for  the  density 
distribution,  N( r,  /),  and  the  dominant  distribution  function,  that 
is  7(r,  s,  t)  along  s  —  so,  the  center  moves  as  vt  So  and  the  B ^ 
in  eq  16  are  proportional  to  t2  at  t  — 1  0.  A  distribution  function 
/(r,  s,  t)  for  s  not  close  to  So  is  small  since  F( s,  so,  t)  ~  U  for 
small  t.  The  center  moves  at  a  certain  direction  with  displace¬ 
ment  proportional  to  vt ,  and  the  B ^  in  eq  16  are  proportional 
to  t2  at  /  — **  0.  These  results  present  a  clear  picture  of  nearly 
ballistic  motion  at  t  — ’  0.  With  increase  of  time,  the  motion  of 
the  center  slows  down,  and  the  diffusion  coefficients  increase 
from  zero.  This  stage  of  particle  migration  is  often  called  a 
“snake- like  mode”. 

With  further  increase  in  time,  the  /th  Legendre  component 
in  eqs  4,  14,  and  17,  exponentially  decays  with  a  rate  related  to 
gi.  The  detailed  decay  rate,  g/,  is  determined  by  the  shape  of 
the  phase  function.  Generally  speaking,  the  very  high  /th 
components  decay  in  a  rate  of  order  of  / is ,  as  long  as  its 
Legendre  coefficient  at  distinctly  smaller  than  21  +  1.  Even  in 
the  case  that  the  phase  function  has  a  very  sharp  forward  peak, 
in  which  there  are  non-zero  a\  for  very  high  /th  rank,  the  a\  are, 
usually,  much  smaller  than  21  +  1.  Therefore,  for  the  distribution 
function  at  time  t  after  the  ballistic  stage  is  over,  a  truncation 
in  summation  of  /  is  available. 

At  large  times,  the  distribution  function  tends  to  become 
isotropic.  From  eqs  19— 21,  the  particle  density,  at  t »  lt!v  and 
r  »  1, ,  tends  towards  the  conventional  diffusion  solution  with 
the  diffusive  coefficient  //3.  Therefore,  our  solution  quantita¬ 
tively  describes  how  particles  migrate  from  nearly  ballistic 
motion  to  diffusive  motion. 
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Figure  1.  shows  the  moving  center  of  a  particle’s  density  function,  Rz 
(eq  20),  and  the  diffusion  coefficients,  Dn  and  Da  (eqs  21),  as  functions 
of  lime,  t. 
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IV.  Discussion 

The  cumulant  expansion  terminating  at  the  second  order  is  a 
standard  method  in  statistical  mechanics,10  which  neglects 
cumulants  higher  than  second  order,  and  leads  to  a  Gaussian 
distribution.  If  we  examine  the  spatial  displacement  after  each 
collision  event  as  an  independent  random  variable,  Ar„  the  total 
displacement  is  XAr,  (/  —  1  ,...N),  with  N  the  number  of  collision 
events,  which  can  be  estimated  by  tlpts.  If  we  define  Y  = 
(A0~1/2£Ar„  the  central  limit  theorem  claims  that  if  N  is  a  large 
number,  then  (Yn)c/(Y2)c  ~  N]~n/ 2,  n  >  3.  Therefore,  the  sum 
of  N  variables  will  have  an  essentially  Gaussian  distribution. 
Therefore,  after  enough  collision  events  happened,  the  distribu¬ 
tions  we  have  calculated  are  accurate  in  detail,  not  just  having 
the  correct  center  and  spread.  At  early  time,  the  particle’s  spread 
is  narrow,  hence,  in  many  applications  the  detailed  shape  is  less 
important  than  the  correct  position  and  correct  narrow  width  of 
the  beam,  because  of  the  finite  resolution  of  detection  devices. 

In  case  a  more  accurate  distribution  at  early  time  is  needed, 
by  use  of  eqs  1 1  and  12  with  its  expression  in  Appendix,  and 
a  standard  cumulant  procedure,  the  exact  higher  (up  to  arbitrary 
mh)  order  cumulants  can  be  analytically  calculated.  Then, 
performing  a  numerical  three-dimensional  Fourier  transform, 
the  particle  distribution  function  accurate  up  to  nth  order 
cumulant  approximation  can  be  obtained. 

In  summary,  we  present  an  analytical  solution  of  the  elastic 
Boltzmann  transport  equation  in  an  infinite  uniform  isotropic 
medium.  Using  a  cumulant  expansion  we  can  analytically 
calculate  cumulants  up  to  an  arbitrary  high  order.  By  terminating 
at  the  second  order,  we  have  derived  an  analytical  solution  of 
the  distribution  function,  eq  13,  and  the  density  distribution,  eq 
19,  with  exact  first  cumulant  (center  of  the  distribution)  and 
exact  second  cumulant  (the  half-width  of  spread  of  the  distribu¬ 
tion).  These  expressions  show  a  clear  picture  of  time  evolution 
of  particle  migration  from  ballistic  to  snake-like,  then  to 
diffusion  regime. 
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Appendix 

In  this  Appendix,  we  derive  an  analytical  expression  of  eq 
12  to  nth  order.  By  defining 

K  =  £|;-xj:ru+,,  -  S[,-x «  =  l>.~n  (Al) 


eq  12  can  be  written  as 

A,...,/*)  =  cxp(~/y)  exp(-glt)F<n>(t)  (A2) 


with 


/»(/)  =  -  />/'"  (A3) 

It  is  easy  to  directly  calculate  eq  A3  for  low  n  orders: 

oh i*  j 
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(A4.1) 
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( b  j  -f  b2)b2b3  (fcj  +  b2  +  b3)(b2  +  /?3)/?3 


(A4.3) 


In  each  step  of  integration,  the  difficulty  is  in  determining  the 
constant  term.  In  the  following  we  prove  that  this  term  is  given 
by  (-1  )nHb„(bn  -f  bn-i)...(bn  4-  bn-\  +  ...  4-  b\)].  Equation  A3 
can  be  written  as 


F(nXt)=  f'Qdt'eb/  F”~  *(/')  (A5) 

Using  integration  by  parts  to  eq  A5,  we  obtain 

F<n)(t)  =  UeK'Fln-'Xt)-  f'dt'  e(b"+h'"'yF<n~zXt')]  (A6) 

bn 


Recursively  applying  eq  A6,  we  obtain 
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bnQ>n  +  bn-i)...(bn  +  bn_i  +  ...  +  bx) 


(A7) 


Equation  A7  provides  formulas  to  recursively  evaluate  eq  12 
up  to  nth  order.  Also,  eq  A7  produces  the  above  mentioned 
constant  term.  An  explicit  expression  of  eq  12  can  then  be 
written  as 


(0  =  exp(-/y)  exp(— g;o£- 


(-if  exp(2X-*+ir) 


m=0 


7=1 


(m) 


(A8) 


with  b„+ 1  =  0,  and 


L/m)  ~  2X’  J-m  or  LT  =  X  bk'J  >  m  (A9) 

k=j  k-m+ 1 
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We  consider  an  analytical  solution  of  the  time-dependent  elastic  Boltzmann  transport  equation  in  an  infinite 
uniform  isotropic  medium  with  an  arbitrary  phase  function.  We  obtain  (1)  the  exact  distribution  in  angle,  (2) 
the  exact  first  and  second  spatial  cumulants  at  any  angle,  and  (3)  an  approximate  combined  distribution  in 
position  and  angle  and  a  spatial  distribution  whose  central  position  and  half-width  of  spread  are  always  exact. 

The  resulting  Gaussian  distribution  has  a  center  that  advances  in  time,  and  an  ellipsoidal  contour  that  grows 
and  chan  lies  shape  providing  a  clear  picture  of  the  time  evolution  ot  the  particle  migration  trom  near  ballistic, 
through  snakelike  and  into  the  final  diffusive  regime. 
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I.  INTRODUCTION 

Scientists  have  tried  for  decades  to  develop  exact  or  ac¬ 
curate  analytical  approximate  solutions  of  the  Boltzmann 
transport  equation  in  various  cases  [1-3].  Any  progress  in 
this  direction  is  a  contribution  to  fundamental  research  in 
non-equilibrium  statistical  dynamics.  An  accurate  analytical 
approximation  may  have  applications  in  a  broad  range  of 
fields,  such  as  the  atmosphere,  medicine,  and  solid  state 
physics.  Photon  migration  in  a  highly  scattering  turbid  me¬ 
dium  is  a  good  example.  The  solution  of  inverse  problems  in 
optical  tomography,  such  as  the  location  of  a  tumor  in  a 
woman’s  breast  from  the  scattering  of  light  pulses,  requires 
the  inversion  of  a  weight  matrix  [4]  obtained  by  convoluting 
two  Green’s  functions  of  the  forward  scattering  problem. 
The  analytical  solution  of  the  photon  diffusive  equation  in  an 
infinite  uniform  medium  has  been  broadly  used  as  a  back¬ 
ground  Green’s  function  [4].  By  introducing  “image 
sources,”  the  solution  can  be  extended  to  semi-infinite, 
slabs,  and  boxes  geometry.  The  diffusion  approximation  fails 
at  early  times  when  the  photon  distribution  is  highly  aniso¬ 
tropic.  Solutions  of  the  diffusion  equation  or  the  telegra¬ 
pher’s  equation  do  not  produce  the  correct  ballistic  limit  of 
light  propagation  [5].  The  Monte  Carlo  method  can  be  used 
to  simulate  photon  migration  at  early  times;  however,  de¬ 
tailed  solution  of  a  five-dimensional  Boltzmann  transport 
equation  using  a  predominately  numerical  approach,  with  the 
resolution  good  enough  to  check  the  analytical  solution, 
leads  to  prohibitive  CPU  times. 

Recently,  Polishchuk  et  al.  [6]  and  Pcrclman  et  ai  [7] 
suggested  different  models  of  photon  migration.  They  used 
the  path  integral  approach  and  the  time-dependent  Green’s 
function  method  to  treat  the  photon  migration  problem.  They 
consider  only  multiple  small-angle  scattering,  based  on  the 
fact  that  the  phase  function  (angular  distribution  of  the  scat¬ 
tering  cross  section)  in  many  media  has  a  very  sharp  forward 
peak.  A  solution  of  the  steady  transport  equation  based  on 
the  small  angle  approximation  was  also  presented  by  Ishi- 
maru  [8].  However,  it  can  be  shown  that  the  transport  mean 
free  path  obtained  by  an  average  of  I  —  cos  0  over  small 
angles  could  be  several  times  larger  than  that  obtained  by  an 
average  over  all  angles.  Thus,  the  small  angle  scattering  ap¬ 


proximation  is  not  quantitatively  correct.  Therefore,  a  proce¬ 
dure  permitting  wide-angle  scattering  is  essential. 

In  this  paper,  we  present  analytical  expressions  for  the 
distribution  function  and  the  density  distribution  of  the  solu¬ 
tion  of  the  elastic  Boltzmann  transport  equation  in  an  infinite 
uniform  medium.  The  phase  function  is  assumed  to  depend 
only  on  the  scattering  angle  />(s,s0)  =  P(s-s0).  Under  this 
assumption,  the  small  angle  approximation  is  avoided,  and 
an  arbitrary  phase  function  can  be  handled.  Our  solution  for 
the  distribution  in  angle  is  exact,  as  are  all  first  and  second 
spatial  cumulants  at  any  angle  as  functions  of  time.  After 
many  scattering  events  have  taken  place,  the  central  limit 
theorem  guarantees  that  the  spatial  Gaussian  distribution  cal¬ 
culated  will  become  accurate  in  detail,  all  cumulants  higher 
than  the  second  approach  small  values  relative  to  the  ap¬ 
proximate  power  of  the  second  cumulant.  At  early  times, 
when  the  errors  would  be  worst,  the  spatial  distribution  func¬ 
tion  at  any  angle  is  quantitatively  accurate  in  the  sense  that  it 
has  the  exact  mean  position  (the  first  cumulant)  and  the  exact 
and  narrow  half-width  of  spread  (the  second  cumulant)  as  a 
function  of  time.  Since  the  inverse  scattering  problem  is 
done  with  instruments  of  finite  resolution,  in  the  presence  of 
noise,  finer  detail  is  lost,  and  the  first  two  cumulants  may 
provide  an  adequate  description  of  the  scattered  beam. 

This  paper  is  organized  as  follows.  Section  II  describes 
the  derivation  of  the  formula,  which  includes  (1)  obtaining 
an  exact  solution  of  the  distribution  in  angle,  (2)  obtaining  an 
exact  formal  solution  in  position  and  angle,  (3)  using  the 
cumulant  approximation  up  to  the  second  order  that  leads  to 
a  Gaussian  spatial  distribution,  (4)  obtaining  exact  first  and 
second  spatial  cumulants  based  on  the  exact  angular  distri¬ 
bution.  Section  III  provides  the  main  results  of  the  distribu¬ 
tion  function  in  position  and  angle,  and  the  density  distribu¬ 
tion  in  position  alone.  Section  IV  makes  a  comparison  of  our 
result  for  the  special  case  of  isotropic  scattering  with  that  of 
the  exact  solution  provided  by  Hauge  [9].  A  discussion  of  the 
effectiveness  of  the  cumulant  approximation  is  presented  in 
Sec.  V. 

IL  DERIVATION 

Without  loss  of  generality,  we  discuss  the  photon  scatter¬ 
ing  problem  with  a  given  light  speed  in  the  medium  c .  Ap- 
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plying  our  result  to  an  another  particle  elastic  scattering 
problem,  with  the  constant  particle  speed  in  the  medium  vis 
straightforward.  The  photon  distribution  function  /(r,s,t)  as 
a  function  of  time  /,  position  r  and  direction  s,  in  an  infinite 
uniform  medium,  from  a  point  pulse  light  source  £(r 
-r0)S(s-s0)S(t-0)  obeys  the  Boltzmann  equation  [3] 

r?/(r,s,/)A?/  +  cs-  Vr/(r,s,f)  +  /xfl/(r,s,r) 

=  fi,j  F(s,s')[/(r,s\/)  — /(r,s,J)]ds' 

+  S(r-r0)S(s-s0)S(t-0),  (1) 

where  fis  is  the  scattering  rate,  jj,a  is  the  absorption  rate,  and 
F(s',s )  is  the  phase  function,  normalized  to  fdsrP(s',s ) 
=  1 .  When  the  phase  function  depends  only  on  the  scattering 
angle  in  an  isotropic  medium,  we  can  expand  the  latter  in 
Legendre  polynomials 


P(s,s')=— 2  a,P,( s  s'). 


(2) 


and  regard  a(  as  known,  either  from  Mie  theory  [10],  or  a 
preliminary  experiment. 

We  first  study  the  dynamics  of  the  photon  distribution  in 
the  light  direction  space  F(s,So,f),  on  a  spherical  surface  for 
s  of  radius  1,  which  is  equivalent  to  the  velocity  space  in  the 
elastic  scattering  case.  The  kinetic  equation  for  F(s,s0,f)  can 
be  obtained  by  integrating  Eq.  (1)  over  the  whole  space  r. 
The  spatial  independence  of  jjls ,  and  F( s,s')  retains 
translation  invariance.  Thus  the  integral  of  Eq.  (1)  obeys 


<?F(s,s0,t)/<?r  +  /^F(s,so,r) 


+  /x, 


F(s, 


, so, f)-J  P( s,s' 


)F(s',s0,t)ds' 


=  S(s-s0)S(t-0). 


(3) 


Since  the  integral  of  the  gradient  term  over  all-space  van¬ 
ishes,  in  contrast  to  Eq.  (1),  if  we  expand  F(s,Sq,/)  in  spheri¬ 
cal  harmonics,  its  components  do  not  couple  with  each  other. 
Therefore,  it  is  easy  to  obtain  the  exact  solution  of  Eq.  (3) 
[11]: 


F(s,s0 ,/)  =  2  — j - exp( -g,t)P,(s- So)exp(  -  /xat), 

l  HIT 

(4) 

where  gj  =  >ctv[  1  ~a//(2/*f  l)].  Two  special  values  of  g/  are 
g0  =  0,  which  follows  from  the  normalization  of  F(s,s')  and 
g  i  =  cl l{ ,  where  /,  is  the  transport  mean  free  path,  defined  by 
It  =  c/[fjLs(  1  —  cos  #)],  where  cos  6  is  the  average  of  s  -  s'  with 
F(s.s')  as  weight.  Equation  (4)  serves  as  the  exact  Green’s 
function  of  light  propagation  in  the  velocity  (or  angular) 
space.  Since  in  an  infinite  uniform  medium  this  function  is 
independent  of  the  source  position  r0,  requirements  for  a 
Green’s  function  arc  satisfied,  especially,  a  Chapman- 
Kolmogorov  condition  is  obeyed:  J<7s'F(s",s',r 
-F)F(s\s,F  —  /0)  =  F(s",s./-r0).  In  fact,  in  an  infinite 
uniform  medium,  this  propagator  determines  all  behavior  of 
light  migration,  including  its  spatial  distribution,  because  dis¬ 


placement  is  an  integration  of  velocity  over  time.  The  photon 
distribution  function  /(r,s,f),  for  the  initial  source  direction 
s0  and  the  source  position  r0  =  0,  is  given  by 


/( r,s,f)  = 


r—  c  j  s(tf)dt ' 
Jo 


S[s(f)-s]}9  (5) 


where  the  angle  brackets  denote  the  ensemble  average  in  the 
velocity  space.  The  first  8  function  insures  that  the  displace¬ 
ment,  r-0,  is  given  by  the  path  integral.  The  second  <5  func¬ 
tion  assures  the  correct  final  value  of  direction.  Equation  (5) 
is  a  formally  exact  solution,  but  can  not  be  evaluated  di¬ 
rectly.  We,  hence,  make  a  Fourier  transform  for  the  first  8 
function  in  Eq.  (5)  and  make  a  cumulant  expansion  to  the 
second  order  [12].  For  an  arbitrary  random  variable, 


(^4)^exp((A))exp((A2)c/2),  (6) 


where  index  c  denotes  cumulant:  (A2)c  =  (A2)*~(A)(A).  An 
exact  result  is  valid  only  if  A  is  Gaussian.  In  the  following 
(B)c  is  called  the  cumulant  of  J5,  while  (B)  is  called  the 
moment  of  B.  Substituting  this  approximation  into  the  Fou¬ 
rier  transform  of  Eq.  (5),  we  have 

/(r,s,0  =  F(s,so,<)^27r^  J  dk 

dt'sa(t')j  | 

-  |  J  V  j'odt"nsa(t’)sfiu")] 

-(/>'«.(, ',}(/>Vo)}).  (7) 

where  T  denotes  time-ordered  multiplication  [13].  Integra¬ 
tion  over  k  in  Eq.  (7)  directly  leads  to  a  Gaussian  spatial 
distribution  displayed  in  Eq.  (10)  below.  Using  a  standard 
time-dependent  Green’s  function  approach,  the  ensemble  av¬ 
erage  of  the  cumulants  in  Eq.  (7)  can  be  calculated.  The 
components  of  the  first  cumulant,  which  is  the  average  center 
position  of  the  distribution,  conditioned  on  s=s0  at  r=0  are 
given  by 


1 

F(s,s0,0 


ds*  F(s,sf 


XsraF(  s\s0,f'). 


(8) 


The  denominator  appears  because  this  is  a  conditional  aver¬ 
age.  The  components  of  the  second  moment,  which  is  related 
to  the  second  cumulant  (average  half-width  of  spread)  of  the 
distribution,  conditioned  on  s  =  s0  at  /  =  0  are  given  by 
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dt"T[s„(t9)sp{t")) 


p /(«!  *S2)—  2 

m 


(l  +  m)\ 


/>/<m)(cos6»,)Pim) 


1 


F(s,s0,t) 


X 


J  ds"F(s,s',t-t')s'aF(s',s",t'-t") 


X(cos02)cos[m(</>|-<^2)],  (13) 

where  %=  1  and  Tjm  =  2(m>0),  P)'"\cosff)  is  the  associ¬ 
ated  Legendre  function.  The  recurrence  relations  of  the 
spherical  harmonics  is  given  by 


X  s"pF(  s",  Sq  ,  t")  +  (t.c.) 


(9) 


cos  9' P\m) (cos  ^')=  27TT[</-m+  1)/>/+'(cos  ^') 


where  (t.c.)  means  the  second  term  is  obtained  by  exchang¬ 
ing  the  index  a  and  /3  in  the  first  term.  Equation  (7)  is  the 
only  approximate  formula  used  in  our  derivation.  Formula 
for  calculating  the  first  two  moments,  Eqs.  (8)  and  (9),  are 
exact.  In  Eqs.  (8)  and  (9),  F(s2,S\  j)  is  given  by  Eq.  (4). 
Since  Eq.  (4)  is  exact,  Eqs.  (8)  and  (9)  provide  the  exact  first 
and  second  moments.  Integrations  in  Eqs.  (8)  and  (9)  are 
tedious,  but  straightforward. 


+  (l+m)P\ ln),(cos0')].  (14a) 

sin  d' P\m\cos  n=^[/><-W) 

<?')]•  (14b) 

The  orthogonality  relation  of  the  spherical  harmonics  is 


ffl.  RESULTS 

In  the  following,  we  set  s0  along  the  z  direction  and  de¬ 
note  s  as  ( 6 \  <f>).  Our  cumulant  approximation  to  the  photon 
distribution  function  is  given  by 


J^d  cos  G'P\m\ cos  O' )P(l?\ cos  6') 

2  ( l  +  m)\ 

=  21+ 1  (/-m)!  S"' ' 


(15) 


/(r,s,f) 


F(s,s0,t) 

(4ir)3/2 


(detS)l/2  CXp 


-4(Zr‘ 


)  a/3 


x{r-rc)a(r-rc)p 


(10) 


Using  Eqs.  (13) -(15)  and  making  integrations,  first  over  <f>\ 
then  over  0\  and  last  over  t\  Eq.  (11a)  is  obtained.  Using  a 
siipilar  procedure,  all  results  in  this  section  were  obtained. 

The  square  of  the  average  spread  width  (the  second  cu¬ 
mulant)  is  determined  by 


with  the  center  of  the  packet  (the  first  cumulant),  denoted  by 
rc,  located  at 


Bap=cGAap-  rcarCp/2,  ( 1 6) 


<  =  c2  AtP /(cos  #)[(/+  l)/(g/  —  g/+l)  +  //(g;  — £/_,)], 

(11a) 

*=C2  AiP\'\cos  0)(cos  (f>) 

xlf(gi-gi-i)-f(gi-gi+i)\’  (Hb) 

where  G  =  c  exp(-^/)//r(s,s0,?),A/=(  1/4 7r)exp(“g/r),g/  is 
defined  after  Eq.  (4),  and 

/(g)  =  [exp(gr)-l]/g.  (12) 

rcv  is  obtained  by  replacing  cos  (f>  in  Eq.  (lib)  by  sin  <fi. 

As  an  example,  we  derive  Eq.  (11a)  as  follows: 

r?  =  — - r  I  dt'  \  ds'  F(s.sr  J  —  t,)s'.F(s\s0,t' ), 

‘  f(s,so,t)Jo  J 


where  F(s2,s{J)  is  given  by  Eq.  (4).  We  denote  s 
=  [sx  ,sy  J  =  [sin  0cos  <f>.s\n  0sin  os  6\  The  spherical 
harmonics  addition  theorem  is  given  by  [14] 


with 


A.-=2  A'P/icosG) 


/(/-l)  ^n)  i  (/+  l)(/  +  2) 
21—  1  E‘  +  2/  +  3 


E?) 


l 2 

+  _ — eO) 


£')J'+  E<4) 


21—1  1  21  +  3  1 


(17a) 


^Vn.vv  A \P /(COS  6) 


i  2' 

(/+ 1  )(/+2) 


l(l~  D 
21-1 


'(i) 


-£i2)  + 


21+3 

(/+1H/+2) 


/(/-l) 
21—  1  1 


-(3) 


21  +  3 
Xcos(2  </>) 

1 


~e\4) 


±2  yA,/f>(cos0) 


E{,1)+- 


21—  1  1  2l  +  3‘ 


7(2) 


E?]-- 


1 


7(4) 


(17b) 


21—  1  1  21+3  ‘ 
where  (+)  corresponds  to  A,,  and  (-)  corresponds  to  Avv, 


3874 


W.  CAI,  M.  LAX,  AND  R  R  ALFANO 


PRE  61 


Axy  =  AVJr=X  — A/P/2)(cos  <9)sin(2<£) 


21 -  1 


E) 


( i) 


+ 


-EJ 


(2) _  . 


1 


?o)_ . 
- 1 


1 


2/  +  3  2/—  1  2/  +  3 

1 


e. 


’(4) 


,  (17c) 


A„=A„  =  X  -A,P}n(cos0)(cos^) 

1 


7  2 

2(1  +  2) 


2U-V>e{1) 


21—  1  ' 


2/+  3 


(2)4-  — 1— P<3)+ - -£i4) 


®)2,+  2r7*',+ 


2/  +  3 


(17d) 


Av;.  is  obtained  by  replacing  cos  <f>  in  Eq.  (17d)  by  sin  <f>.  In 
Eqs.  (17a)-(17d) 

E\l)=[f(g,-g,-2)-f(gi-8i-\)V(8i-i-8i-2)’ 

(18a) 


(18b) 

e{3)= [/(«/— «7- i)-tV(gi-8i- 1 )’  (18c) 

£{4)=[/(^(-^+i)_f]/(«(_«'+'>-  (18d) 


A  cumulant  approximate  expression  for  the  pho¬ 
ton  density  distribution  is  obtained  from  N(rJ) 
=  (S[r—cf'0s(t')dt']),  where  an  average  over  the  angular 
distribution  is  required.  Using  /ds£(s,s',r)  =  exp(— fiat),  we 
have  a  Gaussian  shape 


N(  r,f)  = 


1 


(47rDzzct)'m'4'rrDxxct 


-exp 


(z-/?z)2 


4  Dzzct 


Xexp 


(x2+y2) 


exp  (-fiat). 


(19) 


4  Dxxct 

with  a  moving  center  located  at 

/?z=c[l-exp(-g!7)]/gi 
and  the  corresponding  diffusion  coefficients  are  given  by 


(20) 


D--h 


t 

g\ 


3g|~g2 

HJigi-gi) 


[1  -exp(-git)] 


+  — - r[l-exp(-g20] 

g2(gl-g2) 

—  r— j[  1  _exp(  — git)]2! ,  (21a) 

2g  i  J 


D„=Dvv=L  —  +  -+ - - -[1  — exp(— g|/)] 

3/ 1  g  i  g y( g i  g2> 


i 


g2(gl"g2) 


[I  -CXp(-g2f)]  • 


(21b) 


In  contrast  to  Eqs.  (11)  and  (17),  these  results  are  inde¬ 
pendent  of  g,  for  />  2.  Figure  1  shows  the  moving  center  of 
photons,  R.  [Eq.  (20)].  and  the  diffusion  coefficients,  D:z 
and  Dxx  [Eqs.  (21)],  as  function  of  time,  where  g,  are  calcu- 


F1G.  1 .  The  moving  center  of  photon  density  function  R,  (Eq. 
(20)]  and  the  diffusion  coefficients  D;:  and  Dxx  [Eqs.  (21)],  as  a 
function  of  time  t. 


lated  by  Mie  theory  [10]  assuming  (for  this  figure)  w*ter 
droplets  with  rl\  =  1  are  uniformly  distributed  in  air,  with  r 
the  radius  of  the  droplet,  X  the  wavelength  of  light,  and  the 
index  of  refraction  m=  1.33. 

Each  distribution  in  Eq.  (10)  and  Eq.  (19)  describes  a 
photon  “cloud”  anisotropically  spreading  from  a  moving 
center,  with  time-dependent  diffusion  coefficients.  At  early 
time  f->0,  f(g)^t  +  0(t2)  in  Eq.  (12),  and 
+  0(P)  for  j=  1, 2,3,4  in  Eqs.  (18).  From  Eqs.  (11),  Eqs. 
(17),  and  Eqs.  (20)  and  (21),  we  see  that  for  the  density 
distribution,  N(r,t),  and  the  dominant  distribution  function, 
that  is  /(r,s,f)  along  s=Sq  ,  the  center  moves  as  ct  s0  and  the 
Bap  in  Eq.  (16)  are  proportional  to  t 3  at  r— *0.  A  distribution 
function  /(r,s,f)  along  s^Sq  is  small  since  F(s,Sq,0~' 
when  0.  Its  center  moves  at  a  certain  direction  with  dis¬ 
placement  proportional  to  ct ,  and  the  Bap  in  Eq.  (16)  are 
proportional  to  t 2  at  f— +0.  These  results  present  a  clear  pic¬ 
ture  of  nearly  ballistic  motion  at  0.  Roughly  speaking, 
this  near  ballistic  motion  maintains  its  speed  up  to  Rz 
~0.6 1{  [see  Eq.  (20)].  This  closely  agrees  with  experimental 
results  of  optical  coherent  tomography  (OCT)  [15]  that  the 
range  of  good  resolution  extends  to  about  600  /zm,  in  a  tissue 
of  l~  1  mm.  With  increase  of  time,  the  motion  of  the  center 
slows  down,  and  the  diffusion  coefficients  increase  from 
zero.  This  stage  of  photon  migration  is  often  called  a 
“snakelike  mode.” 

With  further  increase  in  time,  the  Ith  Legendre  component 
in  Eqs.  (4),  (11),  and  (17),  exponentially  decay  with  a  rate 
related  to  g(.  The  detailed  decay  rate,  g/,  is  determined  by 
the  shape  of  the  phase  function.  Generally  speaking,  the  very 
high  Ith  components  decays  in  a  rate  of  order  of  fxs ,  as  long 
as  its  Legendre  coefficient  af  distinctly  smaller  than  2/+1. 
Even  in  the  case  that  the  phase  function  has  a  very  sharp 
forward  peak,  in  which  there  are  nonzero  a{  for  very  high  / th 
rank,  the  a}  are,  usually,  much  smaller  than  21+  1.  There¬ 
fore,  for  the  distribution  function  at  time  t  after  the  ballistic 
stage  is  over,  a  truncation  in  the  summation  over  /  is  avail¬ 
able. 

At  large  times,  the  distribution  function  tends  to  become 
isotropic.  From  Eqs.  (19)-(21),  the  photon  density,  at  t 
>ltfc  and  r>l, ,  tends  towards  the  conventional  diffusion 
solution  with  the  diffusive  coefficient  /;/3.  Therefore,  our 
solution  quantitatively  describes  how  the  photon  migrates 
from  nearly  ballistic  motion  to  diffusive  motion. 
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IV.  COMPARISON  WITH  AN  EXACT  SOLUTION 
IN  THE  ISOTROPIC  SCATTERING  CASE 

A  check  of  our  angular  distribution,  Eq.  (4),  the  first  mo¬ 
ments,  Eq.  (11),  and  the  second  moments,  Eq.  (17),  for  a 
special  case  of  isotropic  scattering  is  performed  by  compar¬ 
ing  with  the  exact  solution  given  by  Hauge  [9]  and  agree¬ 
ment  is  verified.  Hauge  has  provided  an  exact  solution  for 
isotropic  scattering  in  the  form  of  a  Fourier  transform  in 
space  and  Laplace  transform  in  lime,  which  is  given  by 

/kf( s)-  dre~'k  rl(r,s,t),  (22) 


weight.  The  results  are  same  as  Eqs.  (27).  Thus,  by  a  slight 
extension  of  Hauge’s  results  we  verify  the  exactness  of  our 
first  moment  in  the  isotropic  scattering  case. 

For  a  check  of  the  second  moment,  we  notice  that  Eqs. 
(18)  are  obtained  from 

|  dt'cxp(at')  I  df  exp(bt") 

Jo  Jo 

[  ( l/b)[(e<a+b)'- l)/(a  +  b)-(c“'-  1  )/a], 
~\(\la)[(ea,  —  1 )/«  —  /],  a=-b. 

(28) 


with 


7k^(s)  — 


£+/u,  +  /k-cs 

1  1 


M  ,  |kk 

l-rj-j-tan  — 
k|c  £+fx 


- 1 


<5(s— s0) 


477  £+fx+ik-csQ  £+ fi+ik  cs0' 


=  1,2,...  . 

Equation  (4)  in  the  isotropic  scattering  case,  reduces  to 
F(s,s0,t)=  —[l-e-n  +  e-^Sis-So).  (24) 

4  77 


Its  Laplace  transform  in  time  is  given  by 


£[F(s,s0,£)]  = 


1 


ft  S(  S“S0) 


4tt  £+/x 


(25) 


In  the  isotropic  scattering  case,  the  limit  as  a  — >0,  or  b— >0,  or 
both  is  needed. 

Equation  (17a),  without  normalization,  in  the  isotropic 
scattering  case  reduce  to 


•  (23) 

#+ cos  6 

1 

4tt 

[P- 

>  paper. 

c2 

t 

t 

4 - — — 

- e-v 

- e  ~  ‘ 

=  fX,  l 

127T 

A6 

1 

a-T-e' 


+  c2-c-'t'<?(s-s0). 


(29) 

This  moment  based  on  our  method  has  a  Laplace  transform, 
given  by 

rrrc,  .2COs2g+COSg  /*  ,  ^  + 

1  “J  477  £(f+/l)3  1277  ^(f+^)3 


U+M-) 


r  £(s-s0). 


(30) 


If  Eq.  (23)  is  evaluated  at  k=0,  that  means  integration  of 
/(r,s,f)  over  r,  the  result  is  the  same  as  Eq.  (25).  Thus  the 
exactness  of  F(s,So,f)  is  verified  for  the  isotropic  scattering 
case. 

The  first  moments,  Eqs.  (11),  without  normalization, 
[without  divided  by  F(s,s0,/)],  for  the  isotropic  scattering 
case,  reduce  by  our  procedure  to 

[  1  +cos  6 ( 1 


(26a) 


j+rc-^s-so) 


f\~c  sin  0  cos  d>  - — 
4  7T 


l-e-* 


-te~^ 


(26b) 


These  coordinates  of  the  center  have  the  Laplace  transforms, 
given  by 


The  corresponding  result  from  Hauge’s  solution  are  obtained 
by  (l/2)d2/d(-ikz)d(-ikz){Eq.(23)}\k=0,  which  implies 
integration  of  (^/*z)/2  with  /(r,s,f)  as  weight  over  space. 
The  same  result  as  Eq.  (30)  is  obtained.  The  similar  proofs 
have  been  performed  for  Kxx,  Kyy,  Kcxz ,  &y_,  and  Kxyy 
verifying  the  exactness  of  our  second  moments.  In  evaluation 
of  the  value  and  the  derivatives  of  B^{  1 
-(^/|k|c)tan“1[|k|c/(^+>a)]}-1  at  k=0,  we  have  S  =  (£ 
+  /x)/(,  Ba= 0,  Baa=2fic2/[3(2((+ /x)],  and  Ba0=  0  if 
a^/3. 

In  the  above  equations  the  term  related  to  ^“^(s-so), 
has  cumulants  rcz  =  ct  and  2 Azz-(r£)2  =  0.  This  spike  repre¬ 
sents  the  unscattered  part  of  the  light,  which  reduces  its  in¬ 
tensity  as  exp(-^.r).  The  scattered  part  of  light  along  the 
directions  of  s^s()  has  the  correct  mean  positions  and 
spreads,  as  has  been  proved. 


£[K]-c 


1  +  cos  6 


S(s-s0) 


477  ^(f+/x)2+(^+At) 


(27a) 


£[^]  — c(sin  fl)(cos  cf> ) 


1  M 
4t r£(f+/x)2' 


(27b) 


Since  moments  can  be  obtained  by  differentiation  of  charac¬ 
teristic  functions,  we  evaluate  d/d(  —  //:„){£</. (23)}|k=0, 
that  means  integration  over  space  of  ra  with  /(r,s,r)  as 


V.  DISCUSSION 

The  decoupling  of  harmonics  is  valid  only  for  the  angular 
distribution,  F(s.s0,r),  because  in  Eq.  (3)  the  term  such  as 
cs- Vr/(r,s,/)  in  Eq.  (1)  disappears.  This  result  is  available 
only  for  an  infinite  uniform  medium,  otherwise  Eq.  (3)  can¬ 
not  be  derived  from  Eq.  (1).  When  the  spatial  related  distri¬ 
bution,  /(r,s,/),  is  calculated,  the  coupling  of  the  different 
harmonics  remains,  and  is  presented  in  Eqs.  (8)  and  (9), 
through  the  recurrence  relation  of  harmonics,  Eq.  (14),  and 
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explicitly  shown  in  Eqs.  (11)  and  (17),  the  results  of  the  first 
two  moments.  Contrasting  with  the  usual  approach  using  an¬ 
gular  moment  expansion  of  Eq.  (1),  our  cumulant  approach 
has  two  remarkable  features:  (a)  since  the  formula  for  calcu¬ 
lating  cumulants,  Eqs.  (8)  and  (9)  (and  possible  extension  to 
higher  order  cumulants),  use  the  standard  Green  s  function 
approach  without  making  approximation  and  the  Green  s 
function,  Eq.  (4),  is  exact,  the  obtained  cumulants,  as  far  as 
the  nth  order  concern,  are  exact,  (b)  The  cumulants  obtained 
appear  as  the  arguments  of  the  exponential  functions  in  Eq. 
(7),  that  implies  that  an  infinite  series  in  the  usual  angular 
moment  expansion  has  been  included.  Therefore,  even 
though  only  derived  by  terminating  at  the  second  order  cu¬ 
mulant,  the  distribution  function  obtained  has  the  exact  cen¬ 
tral  position  and  the  exact  half-width  as  functions  of  time, 
and  thus  leads  to  the  correct  ballistic  limit  at  f— >0  and  cor¬ 
rect  diffusive  limit  at  large  t.  This  result  is  not  achieved  for  a 
general  phase  function  in  any  known  publication. 

The  cumulant  expansion  terminating  at  the  second  order 
is  a  standard  method  in  statistics  [12],  which  neglects  all 
cumulants  higher  than  second  order,  and  leads  to  a  Gaussian 
distribution.  If  we  examine  the  spatial  displacement  after 
each  collision  event  as  an  independent  random  variable,  Arf- , 
the  total  displacement  is  2Arf(/=  1, . . .  ,N).  The  central 
limit  theorem  claims  that  if  N  is  a  large  number,  then  the  sum 
of  N  variables  will  have  an  essentially  Gaussian  distribution. 
Therefore,  after  enough  collision  events  happened,  the  distri¬ 
butions  we  calculated  become  accurate  in  detail,  not  just 
having  the  correct  center  and  spread.  At  early  time,  the  pho- 
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ton  spread  is  narrow,  hence,  in  many  applications  the  de¬ 
tailed  shape  is  less  important  than  the  correct  position  and 
correct  narrow  width  of  the  beam. 

In  case  a  more  accurate  distribution  at  early  time  is 
needed,  the  exact  higher  (than  second)  order  cumulants  can 
be  analytically  calculated,  and  Eq.  (7)  can  be  extended  to 
higher  order.  Analytical  expressions  for  exact  spatial  cumu¬ 
lants  up  to  an  arbitrary  nth  high  order  have  been  derived,  and 
will  be  presented  elsewhere  [16].  However,  a  closed  analyti¬ 
cal  form  in  space  is  unlikely  to  result,  and  a  numerical  Fou¬ 
rier  transform  over  k  would  be  required.  We  have  therefore 
terminated  the  current  calculation  at  second  order  in  this  pa¬ 
per 

In  summary,  we  have  derived  an  analytical  solution  of  the 
distribution  function,  Eq.  (10),  and  the  density  distribution, 
Eq.  (19),  for  the  elastic  Boltzmann  transport  equation  in  an 
infinite  uniform  medium.  This  solution  is  quantitatively  ac¬ 
curate  up  to  the  second  order  cumulant  approximation  and 
shows  a  clear  picture  of  time  evolution  of  particle  migration 
from  ballistic  to  snakelike,  then  to  the  diffusion  regime.  The 
first  two  position  cumulants  at  any  angle  and  the  angular 
distribution  are  completely  exact  as  functions  of  time. 
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parallel  planar  geometry.  An  example  of  successful  reconstruction  of  four  simulated 
hidden  absorptive  objects  up  to  2cm  below  the  surface  inside  a  human  tissue  like 
semi-infinite  turbid  medium  is  provided  at  the  end.  ©  2000  Optical  Society  of 
America 
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1.  Introduction 


Research  on  near-infrared  (NIR)  diffusive  light  for  biomedical  imaging  and  diagnosis  has  advanced 
over  the  past  decade  because  of  its  potential  to  be  a  safe,  non-invasive,  affordable  and  superior 
diagnostics.1-3  To  seek  a  methodology  which  provides  fast  data  acquisition  and  reconstruction  to 
perform  imaging  with  high  resolution  in  real  time,  a  variety  of  techniques  have  been  explored  includ¬ 
ing  time  resolved  picosecond  pulses,  continuous  waves,  and  diffuse  photon  density  waves  (DPDW). 
Most  methods  reconstruct  three-dimensional  optical  property  maps  (OPM)  by  a  matrix  inversion 
or  iterative  techniques,  or  by  three-dimensional  rendering  of  two  dimensional  projection  images.4-8 
The  difficulty  of  inverting  the  whole  three-dimensional  map  at  once  is  usually  time  prohibitive  when 
the  number  of  volume  elements  involved  increases,  while  three-dimensional  rendering  of  two  dimen¬ 
sional  projection  images  needs  extra  depth  information  of  inhomogeneities  inside  turbid  media  to 
behave  well  and  has  other  limitations.7,9 

In  this  article,  we  introduce  the  theory  of  propagation  of  the  spatial  Fourier  component  of 
the  scattered  wave  field  inside  the  turbid  medium,  and  then  develop  a  new  optical  diffuse  imaging 
methodology  based  on  this  theory,  using  the  two-dimensional  Fourier  transform  of  photon  inten¬ 
sity  on  a  plane  to  detect  inhomogeneities  in  a  highly  scattering  turbid  medium  when  illuminated 
by  a  picosecond  (near)  plane  wave  pulse.  In  such  a  spatial  Fourier  space,  the  picture  of  photon 
migration  is  much  simplified,  in  the  sense  that  different  spatial  frequency  components  of  the  OPM 
(2D  Fourier  transform  on  the  xy  plane)  are  decoupled  from  each  other  and  only  depend  on  the 
corresponding  spatial  frequency  component  of  the  photon  intensity  on  the  detector  plane.  Based 
on  this  observation,  we  obtain  a  super-fast  reconstruction  of  3D  OPM  by  matrix  inversions  of  each 
spatial  component  independently.  The  effect  of  noise  is  explicitly  handled  by  controlling  the  set 
of  spatial  frequency  components  and  the  regularization  parameters  used  in  the  matrix  inversion. 
After  a  rigorous  account  of  the  theory  and  a  brief  description  of  the  algorithm,  an  example  of 
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reconstruction  of  four  inhomogeneities  located  up  to  2cm  below  the  surface  of  a  human  tissue  like 
semi-infinite  turbid  medium  using  backscattered  photons  is  presented  at  the  end. 

2.  Theory 

The  photon  density  </>( r,  t)  at  position  r  and  time  t  in  a  highly  scattering  turbid  medium  can  be 
described  by  the  diffusion  equation 

j^(r,  t)  -  cV  •  £>(r)  V0(r,  t)  +  c^ia( r)<j>(r,  t )  =  S{ r,  t)  (1) 

The  absorption  coefficient  and  the  diffusion  coefficient  D  =  1/3 fi's  where  /z'  is  the  reduced 

scattering  coefficient,  may  depend  on  the  position  in  the  medium,  c  is  the  speed  of  light  inside  the 
medium,  and  5  is  the  source  term  describing  the  density  of  photons  generated  per  second. 

For  the  case  of  a  uniform  medium  and  an  incident  source  5(r,  t)  (5  =  0  when  t  <  0),  the 
incident  wave  field  is  <pi{r,t)  =  f  d3r'  dt'S(r',t')G(r,r',t  -  t')  where  G(r,r',t)  is  the  Green’s 
function  for  the  diffusion  equation  in  a  uniform  turbid  medium.  When  some  weak  inhomogeneities 
(objects  such  as  tumors)  are  embedded  in  the  medium,  we  write: 

Ma,obj('')  =  Ma  "1 

Mi,  obj(r)  =  Ms+Vs(r)  (2) 

where  (ia  and  /i(  are  the  constant  absorption  and  the  reduced  scattering  coefficients  of  the  otherwise 
homogeneous  medium,  //a,obj(r)  and  /i's  obj(r)  are  the  absorption  and  the  reduced  scattering  coef¬ 
ficients  of  the  embedded  inhomogeneity  which  are  spatially  dependent.  Plug  Eq.  (2)  into  Eq.  (1), 
and  note  the  diffusion  parameter  of  the  inhomogeneity  L>0bj(r)  =  D  +  6D( r)  =  l/3/i(  —  6/i((r)/3/j,s2, 
we  have 

r\ 

—<f>(r,t)  -  DcV2(f>(r,t)  +  /J.accf>(r,  t)  =  S(r,t)  +  cV  •  SD(r)V(j)(r,t)  -  cfyta(r)^(r,  t)  (3) 
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The  complete  right-hand  side  of  Eq.  (3)  now  acts  as  the  source  term,  of  which,  S(r,t)  contributes 
to  the  unperturbed  wave  field  (f> o  =  fa(r,  t),  and  the  rest  contributes  to  the  scattered  wave  field, 


<j)s{r,t)  =  (f>(r,t)  -  4>0(r,t) 

=  j  d\'  dt'G{Y,r\t  -  t'){cVr'  ■  5D{r')VT'(l){r'  ,t')  -  c8iia{r')<i)(r'  ,t')) 

=  -  J dV  J  dt'G(r,r' ,t  -  t')6pa(r')c<j>(r' ,t') 

+  j  dh'  J*dt'^^Vr'G(v,r',t-t')-V (4) 

after  partial  integration. 

To  the  first  order  in  the  variation  of  optical  absorption  and  reduced  scattering  coefficients,  we 
can  just  replace  far',  t')  in  Eq.  (4)  by  <pu  i.e.,  the  total  wave  field  is  a  superposition  of  the  incident 
wave  field  fa  and  the  singly  scattered  wave  field  <j>s ■  This  is  the  well-known  Born  approximation. 

Now  consider  the  configuration  of  mostly  studied  parallel  planar  geometry  with  its  boundaries 
at  z  =  0  and  z  =  d.  Its  Green  function  is  thus10 


=  AnDct  exP^~  ^4£)ct  ~  HaCt)Gz(ztz',t),  ( t  >  0) 


(5) 


where  p  =  p '  =  (xf:y'):  and  Gz{z,z* ,t)  can  be  easily  obtained  by  an  image  method.  Its  two 

dimensional  Fourier  transform  on  p, 

G{q,z,p',z',t)  =  J  d2pG(p,  z,  p' ,  z' ,  t)  exp(— iq  •  p) 

=  exp(-iq  •  p'  -  Dctq 2  -  pact)Gz(z,  z\  t ) 

=  G{q,z,z',t')exp(-iq-  p')  (6) 


For  simplicity,  let’s  restrict  ourselves  to  the  case  of  a  pure  absorptive  perturbation  ( 6pa  /  0 
and  5p's  —  0)  and  of  an  incident  pulse  S(r,t)  =  S(p)5(z  —  zs)6(t)  at  this  moment.  The  scattered 
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wave  field  on  a  plane  0  <  z  <  d  is  thus 


MP,  t)  =  -  I  dV  j  d2p"  J*  dt'G( r,  r',  t  -  t')6pa(r’)cS(p”)G(r' ,  p",  zs,t').  (7) 

from  Eq.  (4)  after  replacing  </>  by  fa.  Inside  Eq.  (7),  expand  the  source  term  S(p"),  and  Green 
functions  G(r,r',t  -  t')  and  G(r',p",zs,t')  into  integrals  of  their  Fourier  components,  we  find: 

fa(p,z,t)  =  -^2)3  /  /  dz1  J  d2p"  jf  dt'  J  d2qG{q,z,z',t-t')exp(iq- (p- p')) 

x6pa{p',z')c  j  d2q"S(q")  exp(iq"  •  p")  J  d2q'G(q',z',zs,t')exp(iq!  ■  {p'  -  p")) 

=  ~  ^  J  d2c j  J  d2q'  J  d2q "  dt'  J  dz' exp(iq  ■  p)G(q,z,z',t  -  t')S(q")G{q!,z',zs,t') 

x  J  d2p'dpa(p\z')exp(-ip'  ■  (q-q'))  J  d2  p"  exp(ip"  •  (q"  -  q')) 

=  -  ^4J2)2  /  rf2<l  /  dt>  J  dz>  exP(*q  '  P)<5r(cl>  t  ~  ~  q',  *') 

xSfa'JGCq'jz'jZ*,*')  (8) 


where  5(q,  zs),  dpa(q,  z)  are  two-dimensional  Fourier  transforms  of  the  source  on  the  z  =  zs  plane, 
and  of  the  perturbation  of  absorption  coefficient  over  the  2  =  z  plane,  respectively.  And  at  last,  we 
recognize  the  two  dimensional  Fourier  transform  of  the  scattered  wave  field  (f>s(p,z,t)  on  a  plane  2 
for  the  case  of  a  pure  absorptive  perturbation 


<Mq>  z,t)  =  -■ 


J  d2q!dz'dp,a{q-  q',z')S(q!,zs)  j  dt'G{q,  z,z' ,t  -  t')G{q[ ,  z' , . 


In  a  similar  fashion,  for  the  case  of  a  pure  scattering  perturbation  (8pa  =  0  and  5p's  ^  0),  the 
two-dimensional  Fourier  transform  of  the  scattered  wave  field  is: 


4>s(q,z,t)  =  J  d2q!dz'8p's{q-q',zl)S{q',zs) 

x  [  dt'{q-q'G(q,z,z',t-t')G{q',z',zs,t') 

Jo 

dG{ q,  2,  z',  t  -  P)  dG{q',  2',  zs,t') 
dz'  dz' 


(10) 
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The  general  Fourier  scattered  wave  field  is  the  sum  of  Eq.  (9)  and  Eq.  (10).  Denote  the 


convolutions 


wa(q,q!,z,t] 


z')  =  [  dt'G(q,z,z' ,t  -  t')G{q',z',zs,t') 

Jo 

n  fl  J±ldG( q,  z,z',t-  t')  dG( q',  /,  za,  t') 

z)  =  /„ - a? - a? 


(ii) 


which  are  the  weight  functions  involved  in  the  propagation  of  spatial  Fourier  components  of  the 
scattered  wave  field,  we  have 


Mq ,z,t)  =  J  d2q'dz' 8{la(q-  q' ,z')S(q' ,zs)wa(q,q' ,z,t;  z1) 

+  12Jy2  J  d2q'dz'8fi's(q-  q' ,  z')S(q' ,  zs) 

x{q  ■  q'wa{q,q! ,  z,t-,  z')  +  ws{q,q' ,  z,t;  z')}  (12) 

A.  Incident  plane  wave 

For  the  simple  case  when  the  incident  wave  is  a  plane  wave  pulse,  i.e.,  S( r,  t)  —  SS(z  —  zs)8(t),  such 
that  S(q,zs)  =  47r255(q),  Eq,  (12)  simplifies  to: 


<£s(q,  z ,  t)  —  —Sc  J  dz  {8fia( q,  z  )nia(q,  0,  z, t;  z  )  ~  ^^2  ws{<li  0}z,t\z  )} 


(13) 


The  most  salient  feature  of  the  above  result  is  that  different  spatial  frequency  components  of  8jj,a 
and  8jj!s  are  decoupled  from  each  other  and  the  q-component  of  the  optical  parameters  depends 
only  on  the  the  corresponding  spatial  frequency  component  of  the  scattered  wave  field.  Thus  the 
dimension  of  the  inverse  problem  to  be  solved  later  is  greatly  reduced,  as  is  the  computation  time. 

Let’s  approximate  the  integration  over  z'  by  a  summation,  and  fix  z  =  Zd  at  the  detection 
plane  (omit  z  hereafter),  the  Fourier  scattered  wave  field  on  the  detection  plane 


N  /  %  \ 

<Mq,t)  =  ScAz  ^2[—8p,a(q,Zj)wa(q, 0,t\Zj)  +  — 3  ,2  —  ws(q,0,t;zj)] 
j  —  \  P S 


(14) 
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where  Az  is  the  discretized  step  size,  N  is  the  total  number  of  slices  (layers)  in  the  z  direction 
between  the  source  plane  and  detection  plane,  and  Zj  is  the  2-coordinate  of  the  central  position  of 
layer  j. 

Set  q  =  0  in  Eq.  (14), 


N  rj  /Q  ^  \ 

0s(O,t)  =  ScAzY^[-8fLa(0,Zj)wa(0,0,t;zj)  +  — J—ws(0,0 ,t;zj)], 
j= i 


(15) 


the  zero  spatial  frequency  components  Sj. iQ(0,  Zj)  and  5p.'s (0,  Zj)  can  be  readily  solved  without  doing 
a  complete  reconstruction.  Due  to  the  nature  of  Fourier  transform,  they  just  provide  the  profile  of 
the  amount  of  total  perturbation  of  absorption  and  reduced  scattering  coefficients  per  slice,  i.e., 
the  depth  profile  of  the  inhomogeneities. 

The  whole  three-dimensional  absorption  and  reduced  scattering  coefficients  map  is  thus  con¬ 
structed  through  an  inverse  Fourier  transform  from  all  the  q-components  of  5fia  and  5p,'s  at  different 
depths,  each  of  which  is  solved  independently  from  a  series  of  time  resolved  scattered  wave  field  <f>s 
by  Eq.  (14).  A  schematic  diagram  of  the  procedure  of  inversion  is  shown  in  Fig.  (1).  The  maximum 
spatial  frequency  (cutoff  frequency)  of  the  components  used  in  the  inversion  is  determined  through 
a  signal-noise-ratio  analysis  and  the  regularization  parameter  in  the  matrix  inversion  is  obtained 
by  the  robust  L-curve  method.11 

Both  transmission  and/or  backscattering  image  reconstruction  configurations  can  easily  be 
made  using  Eqs.  (13)  and  (14). 


B.  Near  plane  wave 

When  the  incident  wave  is  not  a  perfect  plane  wave  but  nearly  plane,  we  write  S(r,t)  =  (S  + 
AS(p))8(z  —  zs)S(t)  where  S  is  the  mean  value  of  S(r,t)  on  the  z  =  zs  plane,  and 

S{q,zs)  =  4ir2S{5{q)  +  f{q})  (16) 
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where  /( q)  <C  1  and  /(O)  =  0.  Regard  4n2Sf{q)  as  a  perturbation  to  =  47r2S<5(q),  producing 
extra  deviations  in  5jla  and  Sfi's  such  that 


6p,a{q,z)  =  8p,^(q,z)  +  ea{q,z) 

8p,s{q,z)  =  8p,W{q,z)  +£s(q,z)  (17) 


where  5/xi°^(q,  z)  and  (q,  are  the  solution  of  Eq.  (13)  with  the  same  scattered  wave  field  on 
the  detection  plane  z  =  z&  for  a  perfect  plane  wave  source  S^.  Plug  Eqs.  (17)  into  Eq.  (12),  after 
rearranging  the  terms  and  retaining  only  up  to  the  first  order  of  perturbation  /( q),  we  get: 


-Sc  j dz'{ea(q,z')wa(q,0,z,t-,z')  -  ^jj^ws(q,  0,  z,  t\  z')} 

=  Sc  J  d2q'dz'f{q,)wa(q,q!,z,t-,z,){8fi<£)(q-q!,z')- 
-Sc  J  d2c{dz'f{a[)ws(q,q[,z,t]z')8jj,{s>\q-al,z') 


q.q'^0)(q-q  ',*')■ 


a 


(18) 


Please  note,  the  left  hand  side  of  Eq.  (18)  is  of  the  same  form  of  the  right  hand  side  of  Eq.  (13), 
and  its  right  hand  side  is  a  function  of  known  variables  and  can  be  readily  calculated.  So  ea(q,  z1), 
£b( q,z')  can  then  be  solved  in  a  similar  way  as  in  sec.  (2A). 


3.  Simulation 

For  demonstration  purposes,  consider  a  semi-infinite  turbid  medium  (z  <  0)  with  its  surface  at 
z  =  0  [Fig.  (2)],  whose  absorption  coefficient  //„  =  0. 0033mm-1  and  reduced  scattering  coefficient 
/. i's  =  1.0mm-1. 

Four  absorbing  objects  A,  B,  C  and  D,  each  6.25mm  x  6.25mm  x  3mm  and  with  absorption 
coefficient  //a,0bj  =  0.02mm-1  and  reduced  scattering  coefficient  equal  to  that  of  the  background,  are 
placed  at  depth  7.5mm,  7.5mm,  19.5mm  and  19.5mm  below  the  surface,  and  their  xy  coordinates  are 
(-25, -18.75)mm,  (-12.5, -3.1)mm,  (9.4, 15.6)mm  and  (9.4, 6.25)mm,  respectively.  The  medium 
is  illuminated  by  an  incident  Gaussian  pulse  of  a  shape  of  exp(— p2/2a2)  with  a  =  50mm  inside  an 
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aperture  of  radius  50mm,  propagating  along  negative  z-axis  at  time  t  =  0. 

These  parameters  are  potentially  applicable  to  optical  mammography  of  the  compressed-breast- 
toward-chest  using  backscattered  photons.  A  series  of  measurements  (total  15  snapshots  from  300ps 
to  2400ps)  of  an  area  100mm  x  100mm  on  the  surface  plane  z  =  0  are  computed  by  using  the  direct 
calculation  for  the  Gaussian  pulse  in  r  space.  The  simulated  data  are  used  as  input  for  inversion 
after  adding  a  1%,  5%  or  10%  Gaussian  noise. 

In  the  reconstruction  part,  the  near-surface  region  of  the  turbid  medium  of  depth  up  to  3  cm 
is  sliced  into  IV  =  10  layers,  i.e.,  A z  =  0.3cm,  and  objects  A  and  B  are  then  located  on  layer  3,  C 
and  D  on  layer  7.  The  detection  plane  of  an  area  of  10  x  10cm2  is  divided  uniformly  into  a  32  x  32 
grid.  Objects  A,  B,  C  and  D  all  take  2x2  elements  by  this  grid.  The  results  of  reconstruction  are 
shown  below. 

A.  Depth  profile 

The  absorption  depth  profile,  i.e.  the  total  absorption  perturbation  per  layer  /d2p5/ia(p,z)  vs 
depth  z  is  shown  in  Fig.  (3)  with  different  noise  levels  for  cases  (a)  1%  noise,  (b)  5%  noise  and  (c) 
10%  noise.  In  (a),  there  are  one  peak  at  depth  z  =  0.75cm  (layer  3)  where  objects  A  and  B  are 
embedded,  and  another  peak  at  z  =  1.95cm  (layer  7)  where  objects  C  and  D  are  embedded.  The 
width  of  the  first  peak  at  half  height  is  0.34cm,  about  the  thickness  of  one  layer  (0.3cm),  which 
means  the  depth  of  objects  A  and  B  is  resolved  very  well.  The  second  peak  of  objects  C  and  D 
spans  two  and  a  half  layers  with  its  width  of  peak  at  half  height  0.74cm,  but  its  peak  position  is 
still  correct. 

When  the  level  of  noise  increases,  the  peak  values  of  both  peaks  decreases,  and  the  half  width 
increases.  The  effect  on  the  second  peak  at  z  =  1.95cm  is  much  more  significant  than  that  on  the 
first  one  at  z  =  0.75cm. 
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B.  Layer  reconstruction 


The  full  3D  OPM  is  reconstructed.  The  reconstructed  absorption  coefficient  of  the  layers  at  the 
two  peak  positions  are  shown  in  Fig  (4),  Fig.  (5)  and  Fig.  (6)  for  the  three  noise  levels  respectively. 
Fig.  (4)  shows  objects  A  and  B  are  clearly  resolved  as  two  objects  centered  at  their  original  positions 
with  negligible  expansion;  and  objects  C  and  D  at  depth  z  =  1.95cm  are  also  detected  at  the  correct 
central  positions,  but  the  resolved  images  are  expanded  on  the  xy  plane.  With  the  increase  of  noise 
level,  the  shape  of  objects  A  and  B  blurs  from  Figs.  (4a)  to  (5a)  and  (6a),  and  the  blur  is  even 
worse  for  objects  C  and  D  under  the  same  condition  (from  Figs.  (4b)  to  (5b)  and  (6b)). 

The  reconstructed  absorption  parameter  for  objects  A  and  B  is  0.0071mm"1  at  noise  level  1%, 
about  36%  of  the  original  value  0.02mm"1  of  the  absorptive  inhomogeneity.  In  other  words,  the 
objects  appears  larger  in  space  with  a  weakened  absorption  parameter.  As  the  noise  level  increases, 
the  effect  is  accentuated  with  a  further  reduction  in  the  resolved  absorption  parameter. 

4.  Discussion 

In  the  inversion  part  of  the  reconstruction,  we  have  approximated  the  Gaussian  wave  by  a  simple 
plane  wave.  A  refined  solution  could  be  obtained  as  discussed  in  sec.  (2B).  But  it  is  not  worth  the 
effort,  as  our  calculation  shows  the  correction  is  less  than  1%  at  layer  2. 

The  image  reconstruction  method  provided  is  fast.  The  time  it  takes  to  perform  an  inversion 
in  the  above  example  (of  32  x  32  x  10  volume  elements)  is  less  than  half  a  minute  using  a  scripting 
language  Python  on  one  180Mhz  CPU  of  Origin  200  computer  from  Silicon  Graphic  Inc.  This 
algorithm  scales  only  linearly  with  the  number  of  elements  in  the  xy  grid,  so  it  can  be  used  to 
handle  larger  data  sets  in  real  time  with  little  difficulty. 

This  approach  does  not  limit  the  number  or  the  thickness  of  the  inhomogeneities.  It  allows 
multiple  inhomogeneities  and  one  inhomogeneity  may  span  several  layers. 
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With  little  effort,  a  depth  profile  (the  sum  of  the  perturbation  of  optical  parameter  vs  depth) 
of  the  inhomogeneities  inside  a  highly  scattering  turbid  medium  can  be  obtained.  This  information 
may  be  very  useful  in  some  cases.  When  the  inhomogeneity  is  found  to  exist  only  in  one  layer  from 
the  depth  profile,  the  summation  in  Eq.  (14)  no  longer  exists.  A  direct  inverse  Fourier  transform  can 
thus  be  used  to  resolve  the  inhomogeneity  when  it  is  a  sole  absorptive  or  scattering  perturbation. 
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